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We present a calculation of the hindered Ml T(2 S) —> ^(lS)^ decay rate using lattice non- 
relativistic QCD. The calculation includes spin-dependent relativistic corrections to the NRQCD 
action through 0(v 6 ) in the quark’s relative velocity, relativistic corrections to the leading order 
current which mediates the transition through the quark’s magnetic moment, radiative corrections 
to the leading spin-magnetic coupling and for the first time a full error budget. We also use gluon 
field ensembles at multiple lattice spacing values, all of which include u, d, s and c quark vacuum 
polarisation. Our result for the branching fraction is B(T(2S) —> r/bilS)^) = 5.4(1.8) x 10 -4 , which 
agrees with the current experimental value. 
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I. INTRODUCTION 

Quantum Chromodynamics (QCD) has been accepted 
as the theory describing the strong force of nature ever 
since the discovery of the Jj ip . Since then, there has been 
a long history of using the spectrum and decays of heavy 
quarkonia in order to understand QCD, heavy quarko- 
nia being the ideal theoretical testing grounds when us¬ 
ing potential models, and more recently, lattice QCD. 
Heavy quarkonium states below threshold are very nar¬ 
row, and electromagnetic transition rates are therefore 
significant. Comparing the theoretical and experimen¬ 
tal rates for these decays then provides a very clear test 
of our understanding of the internal structure of heavy 
quarkonia. 

A certain class of electromagnetic transitions between 
quarkonium states, known as hindered Ml transitions, re¬ 
quire a spin-flip between different radial excitations and 
are particularly sensitive to small relativistic effects [1] 
which can illuminate the dynamics of the initial and fi¬ 
nal state systems. These hindered Ml transitions still 
remain a challenge from both the experimental and the¬ 
oretical perspective. Within the bottomonium sector, 
such decays include the T(2 S) —> r]b(lS)^ radiative tran¬ 
sition, where BaBar measured B(T(2S) —> 77*,(15)7) = 
3.9(1.5) x 1CT 4 [2] in 2009. 

On the theory side, hindered Ml decays have been nor- 
toriously difficult to pin down from within a potential 
model framework [1] , where systematic errors are hard to 
quantify and branching fractions ranging from 0.05x 1CV 4 
to 15 x 10 -4 are found. The reasons for the difficulty in 
accurately predicting these decays from within a poten¬ 
tial model will be discussed in Section VI. The continuum 
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effective field theory approach called potential NRQCD 
(pNRQCD) has been used to predict radiative bottomo¬ 
nium decays, including Ml transitions. While these cal¬ 
culations have become quite precise for the allowed IS —► 
IS Ml transitions, the results for hindered Ml transitions 
are dominated by theoretical uncertainties and presently 
can only give an order-of-magnitude estimate [3, 4]. 

Lattice NRQCD is a first principles tool that has been 
systematically improved by the HPQCD collaboration 
and can aid in reliably pinning down this difficult to 
predict decay. Using this formalism, one can accurately 
overcome each of the issues arising from within a poten¬ 
tial model framework. Previous exploratory work on this 
decay in a lattice NRQCD framework was done in [5, 6]. 
We make a number of improvements to those studies so 
that an accurate calculation can be done, complete with 
a full error budget. Some of these improvements include 
using one-loop radiative corrections in the NRQCD ac¬ 
tion and we show in Section V that these decays are very 
sensitive to a subset of these radiative corrections. 

This paper is organised as follows. In Section II we 
set up notation and formulae relevant to this decay, and 
in Section III we give details of the computational setup 
including a discussion of states in NRQCD at non-zero 
momentum. In Section IV the different currents medi¬ 
ating this transition in NRQCD are shown and the per¬ 
turbative calculation of the matching coefficient from the 
leading order current to full QCD is performed. Finally, 
analysis of the Y(2S) —t 775(15)7 decay rate with a full 
error budget is given in Section V. We conclude with a 
discussion in Section VI. 


II. DECAY RATES FOR RADIATIVE 
TRANSITIONS 


BaBar has measured the branching fraction of the 
T(2 S) -A 77b(15)7 decay as 3.9(1.5) x 10~ 4 [2], which 
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when combined with the T(2 S) total width 31.98 ± 2.63 
keV [7], gives the decay rate 1.25(49) x 10~ 2 keV. 
The large errors on the branching fraction are due to 
the difficulty in isolating the small r]b(lS) signal from 
other nearby photon lines (xbj(2P, IP) -+ Y(15)7, 
T(35, 25) —> T(15)y) and from the large background 
in the energy spectrum of inclusive decays [8]. 

We want to perform an accurate and reliable theoret¬ 
ical calculation to compare to this experimental result. 
Computation of the theoretical decay rate requires the 
matrix element of the appropriate operator between the 
Y(25) and rjb(lS) states as input. In a Lorentz invariant 
theory, using the fact that the matrix element transforms 
as a vector under parity (and parity invariance of our 
theory), the only possible decomposition of the matrix 
element is 


(Vb{mS)(k) | f (0)|T {nS) (p, e(p, A))) = 

2 V^(<7 2 ) 


my(nS) “b ^ribimS) 


-e^ vp,T p v k p e{p, A) c 


( 1 ) 


where q is the photon momentum, e(p, X) a is the polarisa¬ 
tion vector of the T („s) and p = k + q by momentum con¬ 
servation. Using time reversal invariance, one can show 
that V^miq 2 ) is rea l [9]. As the Y(25) is a 66 bound 
state, this Ml (spin-flip) transition can occur by flipping 
the spin on either the quark or the antiquark. Since this 
is a symmetric process, the form factor resulting from 
coupling the current to the quark or to the anti-quark 
is then identical. In our lattice calculation we only cou¬ 
ple the current to the quark (c.f. Sec. IV) and actually 
compute 0(g 2 )|lat = Vnm (l 2 )/ 2 ■ 

The decay rate can now be written as 


r(T(2S) -> 7,6(15)7) = 
16«QEDe 2 |q 


(?7i. T( 2S) + m Vb{1S )) 2 


l^r(0)U 2 (2) 


where ctqed is the fine structure constant, e q is the 
quark charge in units of e (i.e., —1/3 for 6-quarks) and 
|q| = (^T(2S) “ m %(is))/ 2w T(2S) by energy conserva¬ 
tion, ensuring that the photon is on-shell with q 2 = 0. 
Thus, from the theoretical perspective, the most chal¬ 
lenging part of calculating the decay rate from first prin¬ 
ciples is computing the single unknown dimensionless 
hadronic form factor Vjy b (q 2 = 0), which encodes the 
nonperturbative effects of QCD. This quantity can be 
calculated in lattice QCD, and this study will focus on 
the computation of b 2 j (g 2 = 0)|i at . 

Using the experimental value of the decay rate men¬ 
tioned above, as well as |q| = 609(5) MeV measured from 
experiment [2] and ctqed = 1/137, we infer 

U 2 ^( g 2 =0)|ex P = 0.069(14). (3) 


This form factor can be directly compared to V^ rib ( q 2 = 
0 ) | iat ■ From now on, we will drop the |i at subscript to 
avoid superfluous notation. 


III. COMPUTATIONAL DETAILS 

A. Second Generation Nf =2 + 1 + 1 Gluon 
Ensembles 

Our calculation uses gauge held configurations gener¬ 
ated by the MILC collaboration [ 0]. For the gauge fields, 
they used the tadpole-improved Liischer-Weisz gauge ac¬ 
tion, fully improved to 0(a s a 2 ). This is possible as 
the gluon action has coefficients corrected perturbatively 
through 0(a s ), including pieces proportional to the num¬ 
ber of quark flavours in the sea [11]. These ensembles 
are said to have 2 + 1 + 1 flavours in the sea, the up 
and down quarks (treated as two degenerate light quarks 
with mass mi), the strange quark, and the charm quark. 
The sea quarks are included using the HISQ formulation 
of fermions [12], fully improved to 0(a s a 2 ), removing 
one-loop taste-changing processes and possessing smaller 
discretisation errors compared to the previous staggered 
actions. 

Five ensembles were chosen, spanning three lattice 
spacing and three values of mi/m s , so that any depen¬ 
dence on the lattice spacing and sea quark mass could 
be fit and extrapolated to the physical limit. Details are 
given in Table I. Due to the computational expense, most 
of the ensembles use heavier mi than in the real world; 
however one of the ensembles used in this study (set 4 in 
Table I) has physical ami/am s , enabling our calculations 
to be performed at the physical point and reducing un¬ 
certainties associated with unphysically heavy sea quark 
masses. 

Successive configurations generated within each en¬ 
semble are expected to be correlated. These autocor¬ 
relations in meson correlators were studied in [13] for the 
ensembles in Table I. There we find that the autocorre¬ 
lations for bottomonium correlators are not appreciable 
and that the configurations can be treated as statistically 
independent. The ensembles have been fixed to Coulomb 
gauge to allow non-gauge invariant smearings to be used, 
helping extract precise results for the excited states in our 
calculation (c.f. Sec. HID). 


B. 6-quarks Using NRQCD 

This study focuses purely on bottomonium processes, 
and information on these processes can be computed on 
the lattice using combinations of 6-quark propagators, 
calculated on the gluon ensembles listed in Table I. As 
the 6-quark has a Compton wavelength of about 0.04 fm, 
these lattices cannot resolve relativistic 6-quark formula¬ 
tions, owing to a > 0.08 fm. However, it is well known 
that 6-quarks are very nonrelativistic inside their bound 
states (v 2 « 0.1), and thus, using a nonrelativistic ef¬ 
fective field theory (NRQCD) for bottomonium states is 
very appropriate. Within NRQCD, with expansion pa¬ 
rameter v (the velocity of the quark inside the bound 
state), one writes down a tower of operators to a certain 
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TABLE I: Details of the gauge ensembles used in this study. 
P is the gauge coupling, ar is the lattice spacing determined 
from the T(2S — IS) splitting [13], where the error combines 
statistics, experiment and the dominant NRQCD systematic 
error. am q are the sea quark masses, N s x Nt gives the spatial 
and temporal extent of the lattices in lattice units and n c f g 
is the number of configurations in each ensemble. We use 
16 time sources on each configuration to increase statistics. 
Ensemble 1 is referred to as “very coarse”, 2, 3, and 4 as 
“coarse,” and 5 as “fine”. 


Set P ar (fm) 

ami 

am s 

am c 

N s x N t 

72 c fg 

1 

5.8 0.1474(15) 0.013 

0.065 

0.838 

16 x 48 

1020 

2 

6.0 0.1219(9) 

0.0102 

0.0509 

0.635 

24 x 64 

1052 

3 

6.0 0.1195(10) 0.00507 0.0507 0.628 32 x 64 

1000 

4 

6.0 0.1189(9) 

0.00184 0.0507 0.628 

48 x 64 

1000 

5 

6.3 0.0884(6) 

0.0074 

0.037 

0.440 

32 x 96 

1008 


order in v allowing for a systematic inclusion of ever- 
decreasing relativistic corrections. This effective field 
theory is then discretised as lattice NRQCD [14]. There 
are a number of systematic improvements which need 
to be made in order to produce highly accurate results. 
These will be addressed shortly. 

We use a lattice NRQCD action correct through 0(v 4 ), 
with additional spin-dependent 0(v 6 ) terms 1 and include 
discretisation corrections. This lattice formalism has al¬ 
ready been used successfully to study bottomonium 5, 
P and D wave mass splittings [13, 15], precise hyperfine 
splittings [16, 17], B meson decay constants [18], T and 
T' leptonic widths [19] and B } D meson mass splittings 
[17]. The Hamiltonian evolution equations can be writ¬ 
ten as 


G(x, t+l) = e aH G(x.,t) 

G(x, f src ) = </>(x) (4) 


with 


e 


-aH 


aH 0 

a5H 


A< 2 ) 

2 arrib ’ 

aSH v 4 + aSH v 6 ; 




( 5 ) 


(A( 2 )) 2 i 

a5H vi = -c lg(a ^ )3 + C2 8(aTO&)2 


VE-E-V 


- c 3 


- C 4 


8 (arrib) 2 

1 


V X E — E X V 


2 arrib 


A< 4 ) (At 2 )) 2 

er • B + C 5 — -C 6 : 


24am.b 16 n(amb) 2 


aSH v 6 = —C 7 


- c 8 


1 


8 (anib) 3 
3 


64(amb) 4 


|a (2) ,(t.b| 

{At 2 ),<r.( 


V X E — E x V 


)} 


C9 


8 (amb) : 


-<r • E X E . 


( 6 ) 


The parameter n is used to prevent instabilities at large 
momentum due to the kinetic energy operator. A value 
of n = 4 is chosen for all anib values. A smearing function 
</>(x) is used to improve projection onto a particular state 
in the lattice data. Using an array of smearing functions 
to improve the overlap with the ground state and the first 
excited state will prove crucial to obtaining accurate re¬ 
sults for the T(2S) -+ r]b(lS)y decay. To evaluate the 
propagator, we use random wall sources that are imple¬ 
mented stochastically with 17(1) white noise, significantly 
improving the precision of the S-wave states [13]. 

Here, arrib is the bare b quark mass, V is the symmet¬ 
ric lattice derivative, with V the improved version, and 
A( 2 ), A^ 4 ) are the lattice discretisations of E^D 2 , EiD 4 
respectively. E, B are the improved chromoelectric and 
chromomagnetic fields, details of which can be found in 
[13]. Each of these fields, as well as the covariant deriva¬ 
tives, must be tadpole-improved using the same improve¬ 
ment procedure as in the perturbative calculation of the 
matching coefficients [13, 20] (thus removing unphysical 
tadpole diagrams from using the Lie group element rather 
than the Lie algebra element in the construction of the 
lattice field theory). We take the mean trace of the gluon 
field in Landau gauge, uql = (|Tr U^(x)), as the tadpole 
parameter, calculated in [13, 18]. 

The matching coefficients Ci in the above Hamiltonian 
take into account the high-energy UV modes from QCD 
processes that are not present in NRQCD. Each c,; can 
be expanded perturbatively as Cj = 1 + c- 1 ^ a s + 0(a 2 ) 
and, after tadpole improvement, we expect the radiative 
corrections c- 1 ) to be 0(1). Each c- 1 ^ can then be fixed 
by matching a particular lattice NRQCD formalism 2 to 
full continuum QCD. These corrections have previously 
been computed [13, 20]. Alternatively, particular c^’s can 
be tuned nonperturbatively, which we discuss in Section 
VB 9. 


The quantities relevant to this study are insensitive to the spin- 
independent 0(v 6 ) terms within our precision. 


Changing the NRQCD action can modify the Feynman rules used 
in the computation of in perturbation theory, in general 
changing its value. 
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A high-precision calculation with a reliable error bud¬ 
get will require knowledge of at least the 0(a s ) correc¬ 
tions to the matching coefficients. For example, when 
tuning the quark mass amb fully nonperturbatively in 
NRQCD, one computes the kinetic mass of a hadron 3 
[13]. This kinetic mass depends on the internal kinemat¬ 
ics of the hadron, and hence on the terms ci, C 5 , and Cg 
in the Hamiltonian. Using the one-loop corrected coeffi¬ 
cients to these terms has a small but visible effect on the 
kinetic masses and hence on the value of the tuned amb 
[13]. 

In addition to this, for an 0(v 4 ) NRQCD action with 
C 4 = 1 , the kinetic mass for the r]b is actually found to be 
larger than that of the T [13], opposite to what is seen at 
zero momentum and, more importantly, in experiment. 
The explanation is that the a ■ B term gives rise to the 
hyperfine splitting, and the splitting from this term is 
correctly included in the static mass (the mass at zero 
energy, offset due to removing the mass term from the 
Lagrangian). However, relativistic corrections to er • B 
(the term proportional to C 7 in the Hamiltonian above) 
are needed to correctly feed this splitting into the ki¬ 
netic mass. On a fine lattice, a value of C 4 = 1.18 and 
C 7 = 1.25 was needed to yield a hyperfine splitting using 
kinetic masses which agreed with experiment within er¬ 
rors [16]. In order to remove the sensitivity to the cr ■ B 
term when tuning amb , one does not use the kinetic mass 
of a single state, but the spin-averaged kinetic mass of the 
T and r]b [13, 21]. Including aSH v s terms in the evolution 
equations makes the rib kinetic mass lower than that of 
the T, as they include relativistic corrections to the a ■ B 
term. The spin-averaged kinetic mass gets smaller and 
the bare quark mass gets larger [16]. 

The parameters used in this study are summarised in 
Table II. There, Ci,Cs and cq are the correct values for 
a v 4 NRQCD action [13], but the small changes to these 
coefficients in going to a v 6 NRQCD action have a neg¬ 
ligible effect on the quantities studied here, as shown in 
Figure 7. While the amb values from ensembles 1,2 and 
5 listed in Table II have all been tuned against the spin- 
averaged kinetic mass using the Hamiltonian above [16], 
the amb values from ensembles 3 and 4 were previously 
tuned without the aSH v 6 terms [18]. Ensembles 2,3 and 
4 are all coarse lattices and only differ by having dif¬ 
ferent light quark masses in the sea. Ensemble 2 has a 
correctly tuned amb = 2.73 for the Hamiltonian we use, 
corresponding to mj = 4.418 GeV. It is appropriate to 
tune the amb values on the other coarse lattices to match 
this physical value. Using the lattice spacings listed in 
Table I, we find the amb values on ensemble 3 and 4 
listed in Table II. All these ensembles have essentially 
the same value of the lattice spacing, so the running of 


TABLE II: Parameters used for the valence quarks, amb is 
the bare 6 -quark mass in lattice units, uql is the tadpole 
parameter. The c; are coefficients of terms in the NRQCD 
Hamiltonian (see Eq. 6 ). Details of their calculation can be 
found in [13, 20]. 03 , 07 ,cs and 09 are included at tree-level. 
We also list the values of a s used to determine the one-loop 
corrections in the perturbative matching in Sec. IV A and for 
the error budget in Sec. VD. 


Set 

amb 

UQL 

Cl, Cq 

C2 

C4 

C5 

a a (7r/a) 

1 

3.31 

0.8195 

1.36 

1.29 

1.23 

1.21 

0.275 

2 

2.73 

0.8346 

1.31 

1.02 

1.19 

1.16 

0.255 

3 

2.68 

0.8349 

1.31 

1.02 

1.19 

1.16 

0.255 

4 

2.66 

0.8350 

1.31 

1.02 

1.19 

1.16 

0.255 

5 

1.95 

0.8525 

1.21 

0.68 

1.18 

1.12 

0.225 


TABLE III: The local bilinear operators used in this study. 
Note the iy 5 is needed to make the overlaps real [9]. The sec¬ 
ond column gives the J PC states that these operators create 
at rest in an infinite volume continuum. The third column 
gives the helicity eigenvalues A that these operators create at 
nonzero momentum in an infinite volume continuum which 
is only rotationally invariant, while the J in brackets are the 
states which contribute to that helicity (c.f. Section HIE). 


G r (x) J PC 


J p ) 

0 

o-(< 

- J p =0-,l+,2-,...) 


0+(< 

- J p =0+1-,2+,...) 

'i/>y 1 '!^ 1 


|i|(* 

- J=l,2,3,...) 


the bare mass is a negligible effect. This was observed 
with a 0(v 4 ) Hamiltonian [13]. 

Within NRQCD, the Dirac field T can be written in 
terms of the quark if) and anti-quark \ as = (4’iX) T - 
The propagator is then found to be 

S{x\y) = ( G ^ x ^ 0 \ 

' \ 0 -G x (x\v)J 

where G^{x\y) is the two-spinor component quark prop¬ 
agator and G x {x\y) is the two-spinor component anti¬ 
quark propagator, y 5 hermicity becomes G^(x\y) = 
~ G x{y\x)- As such, we write our interpolating operators 
as in Table III and then use the above decomposition, 
with suitable boundary conditions, to write the correla¬ 
tor in terms of G^(x\y). 


! The static mass (the energy corresponding to zero-spatial mo¬ 
mentum) in lattice NRQCD [1 ] is shifted due to the removal of 
the mass term from the Hamiltonian and so one can only tune 
static mass differences fully nonperturbatively. 


C. Non-Integer Momentum on the Lattice 

Using periodic boundary conditions (PBC) for the 
quark fields forces the momentum components to be 
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Pi = 27 riii/L, where is an integer. The issue with this is 
that processes which occur at a specific momentum, such 
as that needed for an on-shcll photon in the form factor 
V 21 Vb (q 2 = 0); cannot be reached at an integer-valued 
momentum. Here, we use “twisted boundary conditions” 
(0BC) [22, 23] in order to find the matrix element at the 
physical q 2 = 0 point. There are some subtleties with 
using 0BC in our calculation that, to our knowledge, are 
not found in the literature, and we give an explicit ex¬ 
ample of the construction of our twisted correlators in 
Appendix B. As seen there, and confirmed by numerical 
data, the twisted and untwisted correlator data should 
agree (if the same momentum is used) on a configuration 
level if everything is done correctly. 

In our calculations, we choose Pi = Pf = q = 0 and 
only twist a single propagator so that p® = — q 6 — 6. 
The choice of isotropic twist momentum 6 = Xo(l> 1 ; 1) x 
2-k/L that gives q 2 = 0 depends on the specific process 
under study and for the T(2S I ) —> decay y 0 is 

found from (2) as: 


_ L W T(2 S) 

2\/3it 2mx(2 s) 


( 7 ) 


has the form 

C f 2pt(^src5 Hsriki P i = 

F src S(0\x',n src ' : n sn k)T sn kS (x|0) 


^V zx ' p8 Tr 


where S e is the twisted propagator (c.f. Appendix B). We 
use smearing functions (j> src (r ), cj) snk (r) on the anti-quark 
field at the source and sink respectively. We employ 
hydrogen-like wavefunctions which have been successful 
in previous studies of 6-physics: <j>(r) = S r ^, exp(—r/ro), 
(2ro — r) exp(—r/2ro). ro is the smearing radius, and we 
point the reader to [13] for further details on the smear- 
ings 5 . The different smearing combinations used in this 
study give a 3 x 3 matrix of correlators. We do not smear 
the quark fields due to complications on using twisted- 
smeared fields as outlined in Appendix B. 

The two-point correlator in (8) can be spectrally de¬ 
composed as 


^2pt iP'sro ^suki P 5 f-) ^ ^ Q(^artfei tya {jlsrcity^ 

fc=1 


( 9 ) 


yielding |q°| 2 = |0| 2 . We choose an isotropic momentum 
as it has been shown to reduce discretisation errors from 
rotational symmetry breaking [13]. Since static masses 
obtained from correlators at rest are shifted by an ar¬ 
bitrary value in NRQCD, tuning \ 0 from lattice data 
would require a more lengthy computation of the kinetic 
masses. Instead, we use the experimental values of these 
masses [7] to tune \o and check that q 2 = 0 from the 
results. 


D. Energies and Amplitudes from Lattice QCD 

Extracting matrix elements on the lattice requires 
knowledge of the lattice amplitudes and energies corre¬ 
sponding to the states being studied. The lattice quantity 
which most naturally encodes information on these is the 
two-point correlator 

C*2pt {jlsrci 'H j snk j P ? t) = 

71 e _ix ' pS ( 0(n snk ; x, t + t 0 )O t (n src ; 0, t 0 )) (8) 

X 

Here, to is the source time, n src , n sn k are the smearing 
type (discussed below) and p e is the twisted momentum. 
After performing the Wick contractions with the bilinear 
operators listed in Table III, the connected 4 correlator 


where Ek is the (k — l) th energy excitation of the in¬ 
terpolating operator 0(x) used in the construction of 
the correlator and a(n src /n sn k , k ) are the corresponding 
amplitudes, labelled by the smearing used at the source 
or sink. We are only interested in the first few excited 
states, so we do not need to worry about multiparticle 
states or the open 6-threshold. Our two-point correlators 
are propagated for a maximum of t/a = 15 timeslices, 
as after this the locally smeared correlator on a fine lat¬ 
tice is largely saturated by the ground state. In addi¬ 
tion, correlators were calculated with 16 different time 
sources on each configuration in order to increase statis¬ 
tics. To avoid complications due to correlations between 
these time sources, correlators were then averaged over 
all sources on the same configuration. 

We fit the 3x3 matrix of correlators from t/a = 
1 — 15 using a simultaneous multi-exponential Bayesian 
fit [ 4, 25] to the spectral decomposition in (9). Dif¬ 
ferent smearings give rise to different amplitudes and so 
we take priors on them to be 0.1(1.5). The priors on the 
ground state energies are estimated from previous results 
and given a suitably wide width [13]. For the zero mo¬ 
mentum case, prior information tells us that the energy 
splittings E n+ i — E n are of the order 500(250) MeV, while 
for the nonzero momentum case, priors of 480(250) MeV 
are used (due to the inclusion of additional states in the 
correlator, see Sec. HIE). Logarithms of the energy split¬ 
tings are taken in the fit to ensure that the ordering of 
states is preserved, helping the stability of the fit [25]. 


4 Disconnected diagrams for heavy quarkonia are expected to be 
negligible as they are suppressed by the heavy quark mass [9]. 


5 We use the smearing types Z, g , e as described in that reference. 
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FIG. 1: The first three energies extracted from the lattice 

5 

NRQCD correlator data with the operator 0 1 across multi¬ 
ple momenta. Statistical errors only. At nonzero momentum, 
the energy of the first excited state is lower than the energy 
of the first excited state at zero-momentum. This is a con¬ 
sequence of new states being present in the correlator data 
at nonzero momentum, as described in Section IIIE. Thus, 
care must be taken not to misidentify states. aE exp ' repre¬ 
sents the energy of the states according to a nonrelativistic, 
rotational dispersion relation reconstructed using the experi¬ 
mental masses details of which can be found in the text. 


E. Energy Eigenstates in Lattice NRQCD 

Theoretically, particle states living in the Hilbert space 
are classified in terms of invariant quantities within irre¬ 
ducible representations (irreps) of the symmetry group of 


FIG. 2: As in Figure 1 but with the operator O 1 . 

a theory. For our calculation, two groups need to be con¬ 
sidered: the Lorentz group and the continuous rotational 
group in three dimensions. Appendix A reviews the con¬ 
struction of the irreps of both these groups at zero and 
nonzero momentum. 

As is well known, the irreps of the Lorentz group at 
rest are described by \p 2 = m 2 \J PC ,M) 1 where J, M 
are the total and third component of angular momentum 
respectively. P is the parity quantum number and for 
quarkonia C is the charge conjugation. The quantum 
numbers J PC classify all particles seen in experiment to 
date [7]. 

However, the symmetry group of NRQCD is only the 
rotational group. At zero momentum, the states within 
such a theory are also described by |p = 0; J PC , M). At 
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nonzero momentum, the situation is significantly differ¬ 
ent, and the irreps are described by |p ^ 0; A), where A 
is an eigenvalue of the helicity operator A = p- J/E. This 
has important consequences for the energy spectrum ex¬ 
tracted from our lattice calculation (compare the zero 
and nonzero momentum lattice spectrum seen in Figures 
1, 2 ) and therefore needs to be fully understood in order 
to have a reliable computation. 

At rest the bilinear operators that we use in our calcu¬ 
lation, listed in Table III with T = i 7 5 , 7 % overlap onto 
definite J PC = 0 h , 1 energy eigenstates respectively 
in the infinite volume continuum version of our theory 
(which is rotationally invariant) [26]. In Appendix A (as 
in [26]) it is shown that at nonzero momentum, O 7 (p) 
is a helicity operator which creates a definite A = 0“ 
energy eigenstate, but O' 1 '( p) creates an admixture of 
A = 0 + ,±l eigenstates, where these A get contributions 
from J p values as listed in the third column of Table III. 
The ± superscript on the A = 0 represents the eigenvalue 
fj = P(—1) J from the II symmetry (a parity transfor¬ 
mation followed by a rotation to bring the momentum 
direction back to the original direction) [2 ]. 

In the correlator data from using O' 1 (p ^ 0), 
guided by the experimental masses and this analysis, 
the lowest states in the spectrum should be 77 ^ (15') (= 
0~ + ),Xbi( 1 -P)(= 1 ++ ),T7 b (25)(= 0 h ), etc. whereas 
from using O 1 ' (p) the lowest states in the spectrum 
should be T(1S)(= 1 ), h b (lP)(= 1 +”),T( 2 S)(= 

1 ), etc. These are the J p states which we see in our 

lattice spectrum at nonzero momentum. 

The first three states extracted from the spectrum with 
the operator O 1 , O 1 are shown in Figures 1, 2 respec¬ 
tively. On the same plot, the solid lines represent the 
energy of the states according to a nonrelativistic, ro¬ 
tational dispersion relation reconstructed using the ex¬ 
perimental masses, e.g., aI5(|p|) = am™ + |p| 2 /2a?n kln , 
where m km is the kinetic mass which we set equal to the 
experimental mass, and ro slm is the static mass offset due 
to neglecting the mass term in the NRQCD Hamiltonian. 
We find am slnl in the correlator data from the O 1 op¬ 
erator by taking the ground state lattice energy at zero 
momentum and finding the shift in the static mass as the 
difference a A = am^ is ^ — aTO^ 1 S j. We then use this 

value of the shift to find am e jp£ lm = am e jp C — a A, to be 
used in the above dispersion relation. We found the shift 
in the O 1 correlator data in the same way. 


The important point to observe in these figures is that 
at nonzero momentum the energy of the first excited 
state is actually lower than the energy of the first excited 
state at zero-momentum, opposite to what one would ex¬ 
pect from a dispersion relation. The reason is clear: at 
nonzero momentum energy eigenstates have definite he¬ 
licity, not definite J p . Therefore our correlator data gets 
contributions from the J p states listed in Table III. 

We conclude that, as Figures 1 and 2 show, one has 
to be careful in equating the states found in NRQCD at 
nonzero momentum with continuum J PC quantum num¬ 
bers and also in extracting matrix elements involving a 
state inflight. However, here we only extract excited 
states at zero-momentum in order to avoid unnecess- 
sary complications and to obtain high-precision results, 
which can be muddled when extracting excited states in 
flight due to the addition of extra states in the spec¬ 
trum and their small overlap factors as described in Ap¬ 
pendix A. After our analysis, we can then be sure that 
we have extracted the correct matrix element for the 
T(2S') —> 77b(15')7 decay. 

F. Matrix Elements from Lattice QCD 

The simplest quantity which encodes information on a 
meson-to-meson decay matrix element from within lat¬ 
tice QCD is the three-point correlator 

C3£(n. rc , n snk , p e f = -q e ; t, T) = (10) 

^2 e_ * X ‘ pS (°f(. n snk\ x, T)J n { q e ; y, t)0^ (n src \ 0, 0)) 

x >y 

where O r ' 1 . Of are interpolating operators which create 
the initial state with polarisation m and final state re¬ 
spectively, J ra (q e ;y,f) = ?/dr n (q 0 ;y, t)ip is the current 
which induces the transition with n labelling the polar¬ 
isation of the photon, and the twisted momenta are de¬ 
scribed in Sec. IIIC. The three-point correlator is visu¬ 
alised as in Figure 3 where the three points in lattice units 
correspond to: the source point of the initial particle at 
time t 0 (equal to zero in (10)); the position and time 
of the current causing the transition at (y, t); and the 
position and time of the final state at (x, T). After per¬ 
forming Wick contractions on the three-point correlator 
the connected contribution, written in terms of NRQCD 
propagators as discussed in Section IIIB, is 


C^l(n sre ,n snk 


,p 0 = -q e ;f,T) = -^, 


x,y 


Tr 


TTG x (O\x)T f Gl(x\y)T n (q 0 ;y)G^y\O) 


( 11 ) 


where the twisted propagator G e (x\y) is defined in Ap- is unnecessarily expensive as we can use the sequential 
pendix B. Direct computation of the propagator G(x\y) source technique (SST) [9, 27] to yield the desired prop- 






agator, which only requires one further evolution. There 
are two ways to package the G(x\y) propagator in the 
three-point correlator when using the SST. The first is 
called the fixed current method, which requires the in¬ 
sertion time t to be fixed and for propagator 2 in Figure 
3 to be used as a source for propagator 9. However, this 
method does not scale well and is undesirably expensive 
for relativistic quark formalisms. 

The second approach is called the fixed sink method. 
In this approach, one fixes the sink time T and factorises 
(11) as 


C3& (»w 5 ’W: 5 Pf = -q ;t,T) 


= -£< 


,-iy-e 


Tr 


TY l H e \y\0)T n (q d ;y)G^y\0)\ (12) 


with 


H e (y |0) = (y\x)T\G^x\0) 

X 


where we have written H(y |0) in terms of the twisted 
propagator that satisfies periodic boundary conditions 
and used the fact that T f commutes with the expo¬ 
nential as described in Appendix B. We have also used 
the NRQCD 7 5 -hermicity conditions from Sec. IIIB, 
and used G^{x\y) = —G^(y\x) because G(x\y) = 
(ip(x)i/j^ (y)). We can obtain H e (y |0) by using the twisted 
evolution equations with the source e !X ' p rJ:G^(a;|0). 

Clearly, the two methods should give the same corre¬ 
lator data as they only differ in how G(x\y) is packaged. 
We have checked this numerically and found it to be true 
on any given configuration up to machine precision. As 
the fixed sink method is more cost effective, this method 
was used for the calculation. Our program structure can 
be visualised in Figure 3. Propagator 1 is generated with 
a smeared, random wall source at time to and propagated 
to time T where the sink smearing is applied. H 8 (y |0) 
is found by using the source e* x ' p r^Gj,(a:|0) and evolv¬ 
ing backwards in time using the twisted configurations 
to a time 0 < t < T. Propagator 2 is made from the 
same random wall as 1. We then combine propagator 2, 
H e {y |0) and the current as in (12) to obtain the three- 
point correlator. We use the same 16 time sources as in 
the two-point correlator and prior to fitting, all data is 
translated to a common to = 0. 

The three-point correlator (10) can be related to ma¬ 
trix elements of the current by inserting a complete set of 
states [9]. By doing so, and using the rotational param- 
eterisation of the overlaps as described in Appendix A, 
G^" is seen to be anti-symmetric. We average over the 
six nonzero contributions using an isotropic momentum 
as 


s~iV 

°3pt 


1 3 

£ 

1=1 


tin 


s^inm 
*°3pt • 


(13) 



FIG. 3: Setup for the three-point correlator calculation as 
described in Sec. IIIF. Propagator 1 is the anti-quark and 
£(x) is the random noise source as described in the text. 


In addition, inserting the complete set of states also 
leads to the functional form of the fitting function 



= £ ajrisnk, i)V$b*(n src , f)e Eit e £/(T (14) 

ij 


where a(n sn fe,*) and b(n src , f) are amplitudes from 
the two-point fitting function in (9). The two-point 
and three-point correlators can be simultaneously fit 
to (9) and (14) respectively using multi-exponential 
chained [28], marginalised [29] Bayesian fitting. Chained, 
marginalised fitting has been shown to significantly de¬ 
crease the fitting time and produce reliable, precise and 
accurate results if the data is in the limit of high statis¬ 
tics (Gaussianly distributed) [30]. We check that re¬ 
sults are compatible from both with and without chained, 
marginalised fits on a subset of the data. We use a prior 
of 0.1 (0.2) for all Vff and the same priors for the am¬ 
plitudes and energies as in the two-point fits described 
in Sec. Ill D. For each current, we obtain data for fixed 
T = 9,12,15 and the same 3x3 matrix of smearings as 
in the two-point correlators. This allows accurate extrac¬ 
tions of the matrix element as it includes excited state 
contributions. 

The use of a singular value decomposition stabilises the 
fit and is standard practice in the literature [28]. In our 
Bayesian fit, this is performed by setting a tolerance and 
replacing all eigenvalues of the correlation matrix smaller 
than this tolerance times the maximum eigenvalue to this 
value [28]. By doing so, this leads to larger errors in 
the fit results and so is a conservative step. We use a 
tolerance of 10~ 4 . 

The matrix element for the T(2 S) —► 77 ^( 15)7 de¬ 
cay will be proportional to V,^ rlh . By equating the fit¬ 
ting functions to their continuum correlator counterparts 
with conventional relativistic normalisation, parameter- 
ising our overlaps using rotational invariance with the 
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initial particle at rest, we find 


• Aqed 


iT-yl- 


r-p-V) = («) 

where 9 is the twisted momentum described in Sec. Ill C. 
Since the static masses obtained from an NRQCD cal¬ 
culation are shifted, as explained previously, we extract 
V 21 Vb (? 2 ) fr° m V 21 using the same experimental masses 
as in Sec. II. A nonrelativistic dispersion relation was 
used to find E Vb ns), which is appropriate as shown in 
Figure 1. 


IV. Ml RADIATIVE DECAY CURRENTS 

In order to compute the form factor V^ b {q 2 ), we need 
to choose currents which will induce a hindered Ml ra¬ 
diative decay. Within a nonrelativistic framework, it is a 
standard result in the literature [31-33] that the leading 
order contribution to the matrix element is suppressed 
due to the orthogonality of the radial wavefunctions and 
relativistic corrections are necessary. This suppression 
introduces a sensitivity to a range of effects that we must 
test and quantify in order to perform an accurate calcu¬ 
lation. The first of these effects is the fact that next-to- 
leading order current contributions are appreciable and 
we need to include them. 

As we are using NRQCD to simulate the 6-quark, 
choosing the currents from a NRQCD and non¬ 
relativistic quantum electrodynamics (NRQED) effective 
field theory is most appropriate. This effective field 
theory can be found straightforwardly by extending the 
SU{3) Lie algebra of NRQCD to a SU( 3) x t/( 1) Lie 
algebra to produce NRQCD + NRQED [3]. Then, in 
principle, one could discretise the SU(3) x 17(1) theory 
and choose appropriate currents from the resulting op¬ 
erators. However, this introduces complications, e.g. the 
17(1) magnetic field only decouples from the SU{3) chro- 
momagnetic field to leading order in the lattice spacing, 
resulting in lattice artefact currents which are not present 
in the continuum. Calculating such currents would re¬ 
quire more computational resources and make the com¬ 
putation of the matching coefficients more difficult. 

Instead, we are free to choose the currents from the 
continuum NRQCD + NRQED theory and renormalise 
these. It is important to understand the power count¬ 
ing in the NRQCD + NRQED effective field theory in 
order to choose our currents appropriately. Given that 
NRQCD 4- NRQED is a SU( 3) x 17(1) effective field the¬ 
ory, it has two expansion parameters. For NRQCD, we 
have the standard expansion parameter v , where v 2 ~ 0.1 
for bottomonium. The only scale available for the on- 
shell emitted photon is the photon’s energy |g 7 | ~ 0.6 
GeV. Since the photon’s energy is the difference between 
the masses of two heavy S-wave quarkonia, it is of the 
order |(f 7 | ~ mv 2 ~ 0.4 GeV. Thus we can expand our 
effective field theory in terms of v only. 

We summarise the power counting as 


• Bqed , Eqed ~ Wj\ 2 - 

• The standard QCD power counting rules for the 
QCD fields. 

• The knowledge that when a derivative acts on the 
photon field, it gives a factor of |g 7 | and when act¬ 
ing on the quark field a factor of p q ~ mv (as the 
valence quark knows nothing of the photon momen¬ 
tum in the initial quarkonium rest frame). 

Ordering the operators that induce a Ml (spin-flip) 
transition from NRQCD + NRQED, we find (to next-to- 
leading order for our decay and borrowing notation from 

[34]) 


Of = • Bqedi l> 

lull) 

Owi = uwi 1 a ’ Bqed}4> 

oTYljj 

Os = UJs g 7 ^^ cr ' x Eqed — Eqed X D)ip 


Os2 — ^32 


b 

i3eeb 
64 mf 


{D~, <t • ( D x Eqed — Eqed X D)}ip 
Otot = Of + Ow 1 + Os + O S2 ( 16 ) 


Here, iD = iV + gAQ CD T a are all pure QCD covari¬ 
ant derivatives, fields marked QED (QCD) are the QED 
(QCD) fields and oj 7 , are the matching coefficients needed 
to reproduce full QCD+QED from our effective theory. 
Using the power counting rules above, we find Of ~ u , 
Owi ~ V s , Os ~ v 5 and O 52 ~ u 7 . We can then factor 
out the photon and electric charge in order to derive the 
currents Jfc(q 0 ;y,t) which give the decomposition of the 
matrix element in (1). For example, the operator Of 
gives rise to the current 

jp = '0 t ( cr X^q) /c e _ ' q,x ^ . 

2rrib 

We then write all currents as J fc (q e ;y,t) = 
i/7rfc(q e ;y, t)^ so that Ffc(q e ;y ,t) will be what enters 
the three-point correlator as in (12). We use the ter¬ 
minology that the form factor coming from the current 
Jf is called V^ r "’\F = ujfV^i^If, where the tilde im¬ 
plies we have factored off the matching coefficient from 
the form factor in the numerical calculation and this 
should be applied later in the analysis. Similar nota¬ 
tion is used for the other currents and we refer to V^ Vb \i 
as unrenormalised form factors. The final form factor is 

*21 ~ 2 ^ 1 ^ 1*21 li¬ 

lt should be noted that there are other currents (sup¬ 
pressed by v or a s ) that contribute to this decay and 
which might be of interest, notably, the QCD analogues 
of the Ow 1 , Og-, operators arising from choosing the elec¬ 
tric (magnetic) fields in (16) to be gluon fields and the 
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photon coming from the full SU( 3) x t/(l) covariant at loop level in the full QCD + QED theory. These can 
derivative. Other operators are those which only occur be written as 


OwiQGD = —WWlQCDg^i^ {Aqeu • D + D ■ Aq E D,CT ' g Bq C L> IV’ 
OsQCD = USQCDT—^^f 7 • ( Aqed X g Eq C D ~ gE Q CD X Aq E D)'4 > 

° m b 

Ow2 = uw2 D l er • BqedD 1 iI> 

Op'p = • DBqed • D + D ■ Bqedct • D%1>. 


(17) 


When attempting power counting on the QCD operators 
above, it is helpful to draw the Feynman diagram that 
such an operator would produce. Essentially, we need to 
contract the gluon field with another, producing another 
factor of gv 3 at least [ ]. Consequently these operators 

are expected to be of order a s v 8 at most. We confirm 
numerically that the form factors from these QCD oper¬ 
ators are suppressed as expected and they are negligible 
within the errors of our final results. Since ujw 2 ,<jJ p > p 
occur only at loop level they are suppressed by 0(a a ) 
relative to Owi- We will introduce a systematic error for 
neglected currents in the final analysis. 


A. Matching Coefficients for the Currents 

The matching coefficients, Wj, appearing in the oper¬ 
ators in (16) are needed to take into account the high- 
energy UV modes from processes in the full theory but 
not present in our effective field theory. They have the 
expansion 1 + + O(a^). Here we compute the one 

loop correction to the coefficient ujf from the leading or¬ 
der current. Following this, we estimate the errors from 
neglecting corrections that we do not calculate. 

Our calculation of the one-loop coefficient uj e ' > is very 
similar to the computation of the one-loop correction of 
C 4 in [20]. Following that analysis, by matching the cur¬ 
rent from NRQCD + NRQED to continuum QCD+QED, 
we find 


U! 


( 1 ) 

F 


= b 


(!) 

a,QED 


yNR,( 1) 


ytad,( 1) 


7 NR,( 1) _ y NR,{ 1) 
^2 ^a.QED 


(18) 


where b^Q ED = Cf/ 27t is the coefficient of the first order 
correction to the quark’s magnetic moment, computed 
analytically in continuum QCD following standard tech¬ 
niques. As bcr^QED is UV finite, this allows us to di¬ 
rectly equate results obtained on the lattice to those ob¬ 
tained in the continuum, since the difference between the 
schemes for UV regulation is then irrelevant. In the gen¬ 


eral matching procedure the continuum and lattice IR 
divergences cancel in the computation of the radiative 
correction; here, because of the standard Ward Identity, 
the continuum and lattice contributions to wj)' are sep¬ 
arately finite. 

Z^ R , Z 2 R , Z^q ED are the renormalisation factors of 
the bare quark mass, the wavefunction and the current 
Jf from (16). These are calculated in lattice NRQCD. 
We automatically generate the Feynman rules for a spe¬ 
cific NRQCD action (along with the Symanzik-improved 
gluonic action) using the HiPPy package, then compute 
the numerical evaluation of these diagrams using the HP- 
src package [35, 36]. We use the full v 4 NRQCD Hamil¬ 
tonian with spin dependent v 6 pieces as defined in ( 6 ). 
Computation of Zm R ’ < ' 1 ' > and Z is identical to [20]. 
Zrn. R ' <l ^ will get contributions from mean-field correc¬ 
tions which we denote as Zm d ’ ( ' 1 \ We use the Landau 
mean link Uq ; = 0.750 [37]. For the action that we use, 
the tadpole correction is [ 20 ] 


Z: 


tad 


3 \ ( 2 ) 

(am b y) U ° 


(19) 


The NRQCD diagrams contributing to Z^q EE are 
shown in Figure 4. Since we do not actually include 
the QED field in our calculation, there are no tadpole 
factors from this term. Note that Fig. 4(a) is generated 
by the current coming from 'll:'a • Bqed 4 > /2amb being 
inserted at the vertex, and Figs. 4(b), 4(c) and 4(d) arise 
from mixing effects from the higher order currents (that 
we include in the calculation of the decay rate) from 
(16). Computation of the Feynman diagrams shown in 
Figs. 4(b), 4(c) and 4(d) is more involved than that of 
Fig. 4(a), so they are not included in this calculation, but 
we plan on computing them in future work. For now, we 
will introduce a systematic error from neglecting these 
contributions. 

A breakdown of the numerical values of the various 
terms that enter u> E ^ for the masses that we use in this 
calculation is shown in Table IV. tUp 1 was computed for 
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FIG. 4: Classes of one-loop diagrams which contribute to 
Z^qed as described in the text. The cross inside a circle 
represents the Jf current obtained from (16), while the solid 
box represents the higher order currents from (16) and the 
exchange of a gluon is denoted by a curly line. 

TABLE IV: Breakdown of the different terms that go into 
4 1J . The a s (q* = if/a) values are taken from Table II. 


arrib 

1.90 

2.70 

3.30 

7 ( 1 ) 1 7 ( 1 ) 1 7 ( 1 ) 

Zjm -T A 2 ' ^cr,QED 

rytad 

. ,(!) 

UJ p 

a s (7r/a)4 

1.2961(5) 

-1.1233 

0.0394(5) 

0.0089(1) 

0.9061(4) 

-0.8086 

0.1148(4) 

0.0293(1) 

0.7585(6) 

-0.7066 

0.1603(6) 

0.0441(2) 


TABLE V: Values of at various amt values. 


amb 

1.1 

1.5 

2.1 

, , (i) 
LO p 

- 0 . 211 ( 2 ) 

-0.030(1) 0.0626(9) 

amb 

2.4 

4.0 

4.6 

uj p 

0.0918(7) 

0.2039(4) 

0.2372(5) 



a range of masses (neglecting the mixing down) and we 
give these values in Table V. 

We show the values of ujp' 1 with a smooth interpolat¬ 
ing curve in Figure 5. This interpolating curve was cho¬ 
sen to be a polynomial in 1/arrib in order to reproduce 
the static limit as mb —> oo. To fit these values easily 
we increased the errors on the points returned by HPsrc 
to 1%. We use a Bayesian fit to all points in Figure 5 
against a polynomial in 1/amb- We found the smallest 
X 2 /dof(dof) = 0.7(9) and largest Gaussian Bayes Fac¬ 
tor [24] when including all terms in the polynomial up to 
and including the quartic term. We used a prior for the 
constant piece as the polynomial of 0.4(2) and priors for 
the coefficients of the 1 /(amb) n pieces of 0(1). 


B. Systematic Error from Current Matching 
Coefficients 

We need to include a systematic error from not know¬ 
ing the matching coefficients in the currents to infinite 
precision. There are two distinct types of errors in this 
case: the first is from neglecting the 0(a 2 ) corrections 
in wp and the 0(a s ) corrections to the matching coeffi¬ 
cients of the other currents; the second is from neglecting 
the mixing down effects on the values of used in the 
calculation. We will estimate each of these in turn. 

To estimate the effect of neglecting the higher order 
corrections that we have not calculated, it is helpful to 
compare our result to the pure NRQED calculation of 
[34]. There, the authors find that the continuum QED 
contribution to their u>p^ is the anomalous magnetic mo 


FIG. 5: The values ’ with a smooth interpolating curve as 
described in the text. 


ment of the electron a/27r, while for us it is the anoma¬ 
lous magnetic moment of the quark CfOl s / 27t. For the 
NRQED contribution, they find no IR log nor a constant 
piece and in their continuum approach the UV power 
law divergences may be omitted. Although we find no 
IR log in our data, we cannot neglect the UV power law 
divergence associated with the momentum cutoff. This 
shows up as a polynomial in l/amb as mentioned above. 
We observe that this lattice artefact contribution gives a 
negative contribution to the continuum value, as shown 
in Table IV, and for the amb range that we are interested 
in |a s a4| < CfOLsI^I’k. As we are observing similar be¬ 
haviour over this mass range as the pure NRQED calcu¬ 
lation, we can use that calculation to estimate the error 
conservatively. 

As shown in [34] and confirmed by the small values of 
our numerical data, the matching coefficients can actu¬ 
ally be expanded in a s / it. In principle, the second order 
coefficient of ojp could be 0( 1), and then this contribu¬ 
tion could be 0(a 2 /Tr 2 ). As such, we allow for an addi¬ 
tive systematic error (assumed to be correlated across all 
ensembles) of 1 ± a 2 /n 2 from not knowing higher order 
contributions to top- 

We have not included the 0(ot s /7 r) contributions to the 
other matching coefficients in (16), namely uis , u>wi and 
ujs 2 - A difficult calculation would be necessary to de¬ 
termine them. Again, we use the equivalent parameters 
from the pure NRQED calculation [34] to estimate the 
systematic error. The pure NRQED equivalent of %i 
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TABLE VI: Values of the unrenormalised form factors 
V^U, as described in Section V, from the lattice NRQCD 
data on the ensemble labeled set 1 in Table I. We also give ele¬ 
ments of the correlation matrix. A value of a 2 q 2 = 0.0034(21) 
was found from the data. 


p 

Value 

C(p,V |f) 

C(p,V\wi) C(p,V\ s ) 

\/ r Vb I 

V 21 k 

0.1818(42) 



V^ b \wi 

-0.0594(12) 

-0.4010 


V 2 T i Vb \s 

-0.0339(17) 

-0.2932 

0.1261 

V^lsi 

-0.0037(3) 

-0.0624 

0.3488 -0.2678 


TABLE VII: Values and correlation matrix elements of the 
V 21 rlb \i from the ensemble labeled set 2 in Table I. A value of 
a 2 q 2 = 0.00338(92) was found from the data. 


P 

Value 

C(p,V\ F ) C(p,V\wi) C(p,V\ s ) 


0.1765(22) 



F 2 ?Vi 

-0.0593(7) 

-0.5298 


V^ b \s 

-0.0293(8) 

-0.3803 

0.3065 

V 2 T i Vb \si 

-0.0045(2) 

-0.0134 

0.3962 -0.2264 


has log contributions in its first order coefficient and so 
we allow for an additive correlated systematic error of 
1 ± a s /iT to the tree level value. We allow the same error 
on ujs 2 - 

The one loop correction of the pure NRQED equivalent 
of cos was found to be 20 ^!"* = a/ir. As such, we allow 
for an additive correlated systematic error on our cos of 
1 ± CfC( s /w , to compensate for using the tree level value 
in the calculation of the decay rate. This is a conser¬ 
vative estimate as we see above that the lattice artefacts 
actually subtract away some of this contribution over the 
mass range we are interested in. 

The mixing down effects from diagrams (b), (c) and (d) 
in Figure 4 are difficult to estimate since each graph by it¬ 
self can be IR divergent but is IR finite. We allow an 
uncertainty of 30% in the one-loop coefficient (correlated 
across all lattice spacings) from neglecting the mixing 
down. There is no substitute for the actual calculation 
though, and we intend to do this in the future. 

V. RESULTS FOR THE T(25') ->■ ^(lSb DECAY 

The unrenormalised form factors, V^ v (q 2 = 0)|i, for 
each of the currents obtained from (16) are computed for 
each of the ensembles listed in Table I and their values 
are given in Tables VI, VII, VIII, IX and X. A visual rep¬ 
resentation of V^i^q 2 = 0)|i is shown in Figure 6 . From 
this, we can see that the form factor from the current Jf 


TABLE VIII: Values and correlation matrix elements of the 
V^ Vb \i from set 3 in Table I. A value of a 2 q 2 = 0.0007(12) 
was found from the data. 


P 

Value 

C(p,V\ F ) C(p,V\wi) C(p,V\ s ) 

f/ r Vb I 

v 2 i k 
V 2 i lb \wi 
V 2 r i lb \s 
V 2 r r\si 

0.1720(36) 

-0.0577(10) 

-0.0309(12) 

-0.0032(3) 

-0.2634 

-0.1887 0.2733 

0.0213 0.1346 -0.1634 


TABLE IX: Values and correlation matrix elements of the 
V^ Vb \i, from set 4 in Table I. A value of a 2 q 2 = 0.00066(70) 
was found from the data. 


P 

Value 

C( P ,V\f) C(p,V\ wi) C(p,V\ s ) 

f/ r Vb I 

v 2i k 

0.1710(27) 



v^ b \wi 

-0.0596(7) 

-0.4441 


v£">\s 

-0.0289(10) 

-0.3281 

0.2708 

v^ b \si 

-0.0038(2) 

0.0206 

0.1493 -0.3195 


is leading order, and the other currents give a negative 
contribution to Jf of approximately 30%, 20%, 3% for 
Jwu Js , Jsi respectively across all ensembles. Note that 
these values do not appear to obey the power counting 
for the currents given in Sec. IV; however, we understand 
(and explain below) that the leading-order contribution 
is suppressed for these hindered transitions. Similar be¬ 
haviour was seen in previous lattice NRQCD studies of 
this decay [5, 6 ]. 

We also need to determine the sensitivity of our form 
factors to the different parameters used in our calculation 
and use this analysis to give a reliable error budget. This 
is easily done in lattice NRQCD, as we can simply change 
the value of a single parameter and rerun the whole cal¬ 
culation. The results are shown in Figure 7, where we 
denote p as a parameter to vary (either Ci or mb) and use 
A = p test — p°( a s) to signify an upwards/downwards shift 
from the 0(a s ) correct value p°( a ^ (amb is tuned fully 
nonperturbatively but we use amb = p°^ a ^ to avoid ad¬ 
ditional superfluous notation). The values of the changed 
parameters are given in Table XI. 

From Figure 7 we can see that the form factor is most 
sensitive to the value of C 4 , while C 7 and mb are also im¬ 
portant. We need to describe this sensitivity in order to 
give a reliable estimate on the error from not knowing 
each of these parameters to infinite precision. Interest¬ 
ingly, it is useful to note that the sensitivity to these 
parameters comes from the Jf current, as shown in Fig¬ 
ure 8 . We will use a simple potential model analysis to 
understand the deficiencies in the naive power counting, 
where these sensitivities arise from, and to gain insight 
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TABLE X: Values and correlation matrix elements of the 
V^ Vb \i, from set 5 in Table I. A value of a 2 q 2 = —0.0021(6) 
was found from the data. 


p 

Value 

C(p,V |f) C(p,V\ wl ) C(p,V\s) 

V 21 \F 

V^lwi 

V 2 T i Vb \s 

V^lsi 

0.1785(31) 

-0.0618(15) 

-0.0276(10) 

-0.0060(5) 

-0.0703 

-0.0925 0.1526 

0.0457 0.3260 -0.0266 



Jf Jw i Js /si 

Currents 


FIG. 6: The value of the unrenormalised form factor, as de¬ 
scribed in the text, arising from each current across the differ¬ 
ent ensembles listed in Table I. Statistical error only ( « 2—3% 
for each current). 


into this hindered Ml decay. 


A. Phenomenological Insight: Potential Model 
Analysis 

In a potential model framework one would consider pe¬ 
riodic harmonic time-dependent perturbations and find 
the matrix element as the overlap between the spatial 
part of the potential and the initial and final states un¬ 
der study. For an Ml decay, mediated by either of the 
constituent quarks’ magnetic moment er • B, one can find 
the matrix element as [38] (labeling the spatial part of the 
potential as Jf, similar to the current we use in Section 


TABLE XI: Values of the varied parameters used to obtain 
Figure 7. A > 0 (A < 0) denotes an upwards (downwards) 
shift in the parameter as described in the text. for 

A = 0 values are taken from Table II and reproduced here for 
convenience. 


Parameter 

ptest f or A < 0 

for A = 0 

ptest fo r A > 0 

Cl = C6 

1.00 

1.31 

1.50 

C2 

0.75 

1.02 

1.25 

C3 

0.75 

1.00 

1.25 

C4 

1.00 

1.19 

1.50 

C5 

1.00 

1.16 

1.50 

C7 


1.00 

1.50 

m b 

2.5935 

2.73 




FIG. 7: The variation of the form factor with the parameters 
used in this study. A > 0 (A < 0) denotes an upwards 
(downwards) shift in the parameter as described in the text, 
and the values of the varied parameters can be found in Table 
XI. The data for A > 0 (A < 0) were generated on a subset of 
400 configurations of the coarse lattice denoted Set 2 in Table 
I. Statistical errors only. 

IV to highlight comparisons) 

(r] b (mS)\J F \Y(nS)) = 

S f i dr r 2 R* m rib (r)j 0 Rn,r(r) 

with the integral expanded as 

dr r 2 R* m ^ b {r)j 0 (^-^j R n .r{r) = 

Sum + a 2 \q 1 \ 2 rl + a 4 |g 7 | 4 ro H-. (20) 

Here, we have factored the spin piece Sfi in the ma¬ 
trix element from the radial integral (appropriate in the 
nonrelativistic limit) and used the Taylor expansion of 
jo(x ) = sin(a:)/x = J2 n {—l) n x 2n /{2n + 1)! to see that it 
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FIG. 8 : How each of the unrenormalised form factors from the 
different currents vary with C 4 . As can be seen, the sensitivity 
comes from the Jf current. The reason for this is described 
in Section V B 1. 


is a polynomial in |g 7 | 2 . Additionally, the only scale in 
the wavefunctions capable of being combined with |g 7 | 2 
to make it dimensionless is some combination of the Bohr 
radii of each state, which we call ro- The <22/ are coef¬ 
ficients which could be calculated if wave-functions were 
supplied. The leading Kronecker 5-function in (20) comes 
from noting the orthogonality condition in the extreme 
nonrelativistic limit, |g 7 | 2 —> 0. 

As can be seen, for a nS —¥ nS transition, the leading 
order term in (20) is one. However, for transitions be¬ 
tween different radial excitations, the S nm vanishes and 
we are left with a leading order |g 7 | 2 rg term. The radii 
of the bottonronium states under study are of the order 
the reciprocal of the typical momentum, e.g, tq ~ 1 jmv. 
Thus, as |g 7 | 2 7g ~ m 2 v 4 / (m 2 v 2 ) ~ v 2 , the leading or¬ 
der matrix element from Jp in a radially excited decay is 
suppressed by a factor of v 2 more than naively expected 
from using power-counting rules on the currents alone. 
This suppression leads to an array of sensitivities that 
make this decay particularly difficult to pin down the¬ 
oretically from within a potential model [1], as we will 
expand upon in Section VI. 

Due to the derivatives in the other currents listed in 
(16), the matrix elements of these currents give rise to 
wavefunction overlaps that are not orthogonal in the ex¬ 
treme nonrelativistic limit, and as such are not more sup¬ 
pressed for radially excited transitions. The derivatives 
act on the initial bottomonium state and give a lead¬ 


ing order p ~ 0(mv ) effect, which does not depend on 
the photon momentum, as can be seen by taking the 
|g 7 | —> 0 limit. This results in the relativistic corrections 
to the leading order Jp current, which we have included 
in our calculation, having appreciable effects (see Fig. 6), 
namely Jw 1, Js- The orthogonality of the radial wave- 
function muddles up the power counting of the first few 
currents, but additional derivatives in relativistic correc¬ 
tions to these currents will suppress them further. By 
including the current Js 2 , we check that added derivi- 
tives do suppress the contribution of the current further 
as expected. 

By examining (20), we found that the leading order 
matrix element for the radially-excited radiative transi¬ 
tion can be suppressed more than we would naively ex¬ 
pect from just power-counting the current alone. Rela¬ 
tivistic corrections to the Jf current are then apprecia¬ 
ble, explaining the behaviour seen in Figure 6. Even if 
we included the relativistic corrections to the current in a 
potential model, we still would not get the correct value 
for this decay, as we also need to consider all relativistic 
corrections to the wavefunctions arising from perturba¬ 
tive potentials in the Hamiltonian. This gives rise to the 
sensitivities to the different parameters as seen in Figure 
8, which we explain below. To do so, it is sufficient to 
consider first order time-independent perturbation the¬ 
ory. 

B. Sensitivity and Errors from Terms in the 
NRQCD Action 

We want to consider potentials arising from relativistic 
corrections in the NRQCD action causing perturbations 
of the wavefunction. To first order in a s we have 

V Vb 

1*76(1 S)}« = MlS))<°> - £ \ Vb (mS))M-^- 

m^l ml 

T/T 

|T(2S'))( 1 ) = |T(2S')) (0) - ^2 |T (nS)) {0) ^ . (21) 

n^2 En2 

The state In/ 1 ) ( | n))°)) is the first-order perturbed state 
(the unperturbed state), V nm = (°)(n|Vjm/°) with V 
being the potential representing the perturbation and 
E nm = En — Em\ Now, we take currents of interest 
between these states to yield 

W(7 ?b (lA)|J l |T(25))W = 

(°)(7 7b (15)|J l |T(25)) (0 ) 

_ T /Vb * 

-E#- (0) (ilbimS) | Ji |T(2S’)) 

m/l ml 

-Y,3 l0 \*mjmns)r. (22) 

n^ 2 

As mentioned above, for the current Jp , due to the 
fact that (°)(77i,(lS)| Jf\T(2S))^ is suppressed for radi- 
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FIG. 9: The C4 dependence of V^ r,b \f as described in the 
text, along with the lattice values of V^ Vb \F- 

ally excited decays, the ^(r]b(nS)\Jp\T(nS)) < ' 0 ' > pieces in 
the second term in (22) become appreciable. The matrix 
elements arising from currents with derivatives are al¬ 
ready suppressed, and the first order corrections to these 
matrix elements are not appreciable, as seen in Figure 8. 


1. Sensitivity and Error from c 4 : 

Including a potential from the exchange of a single 
gluon between two vertices involving the chromomagnetic 
operator as shown in Appendix C, we find 

«(r ?6 (15)|J F |T(25))W = 

^( Vb (lS)\J F \T(2S))^ + -^-^( 0 )^( 0 ) 

9771^21 

x (6 ( 0) ( % (25 )|J f |T(25))(°) 

+ 2 (°)(?7 6 (15)|J f |T(15))W +0(v 2 )) 

= <°><»fr(lS)|J F |T(2S)>< 0 > 

+ (0)fc(0) + 0(v 2 ) . (23) 

The reason for the sensitivity to c 4 is clear. The matrix 
element ^°^(77f ) (15')| J F |T(25 , ))^°^ is suppressed due to the 
orthogonality of the radial wavefunctions in (20), while 

(i7&(n5)| Jp|T(nS))(°) is not. This results in the second 
term in (23) being sizeable compared to the first. 

Since we have values of the form factor at three values 
of c 4 on a coarse lattice as shown in Figure 8, and an un¬ 
derstanding that the functional dependence of the form 
factor on c 4 should be Vp = a Ci + c 2 b C4 from (23), we 
should check that this is consistent. We use the c 4 = 1.00 
and c 4 = 1.19 values from our lattice NRQCD calcula¬ 
tion listed in Table XII to find the values of a r , and b r ., 
in Table XIII. 

We can also relate the second term from the leading or¬ 
der approximation in (23) to quantities that are measured 


in experiment and check the consistency of the value of 
b Ci given in Table XIII. By comparing the decay rate for¬ 
mulae from a potential model calculation [9] with the one 
given in (2), we find: 


v r Vb _ ( m r(2S) + m ri b (is) \ 

21 V 2 m b ) 


X Jo ^ r2R * m ^ r ^° (^ 2 ~) 

and then using this in (23) yields: 

c\b a = ( "T, 2 S )+m„ ( , s , \ vggggM + 0[v 2 ) 


2m b 


(24) 


where A (iS) is the hyperhne splitting between f’th radial 
excitations. Using the values of c 4 , a and arrib from set 
2 in Table II, along with the PDG average [7] values for 
A (iS) and the spin averaged U 2 1, we find b C4 = 0.105(14). 
This is consistent with the value of b C4 from Table XIII. 

In Figure 9, we show the strong c 4 dependence of 
V**\ F = a Ci + c 2 b Ci , along with the the lattice values 
of V£*\ F shown in Figure 8. This illustrates both the 
need for at least the C2(a s )-correct value of c 4 and the 
consistency of a Ci and b C4 with all our lattice data. 

Since we only know c 4 to one loop in perturbation the¬ 
ory, there will be a systematic error associated with not 
knowing it to higher orders. With the above functional 
dependence of V^ rib \ f = a C4 + c 2 b Ci , an error of 2 a 2 b C4 
should be introduced from not knowing c 4 to second or¬ 
der. As there is little lattice spacing dependence in the 
unrenormalised form factors as shown in Figure 6, we use 
the value of b C4 from Table XII across all ensembles and 
introduce an additive systematic error (correlated across 
lattice spacings) of 2 a 2 b Ci from not knowing c 4 to more 
than one loop. We also allow for the statistical error in 
the determination of c^ coming from the Vegas integra¬ 
tion [20] by adding an error of 2a s 8c^b Ci . 

With the other currents that have derivatives, the sit¬ 
uation is significantly different. Due to the derivatives, 
the second term in (22) is always suppressed and rela¬ 
tivistic corrections are not an appreciable effect, as seen 
in Figure 8. Variations of these currents with c 4 are not 
appreciable within the other errors. 


2. Sensitivity and Error from cj: 

The C'j operator is a D 2 correction to the c 4 term and 
is expected to be a 0(v 2 ) effect. We can proceed as 
before, assuming a linear functional dependence on C7 as 
v£> b \ F = a C7 + C7 b Cj , coming from the exchange of a 
single gluon from a c 4 vertex and a C7 vertex. Using our 
data points in Table XII, we find a C7 , b C7 listed in Table 
XIII. 

It is seen that b C7 gives a negative contribution as a 
consequence of the D 2 and the ratio b C7 /b Ci = —0.20(18) 
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TABLE XII: Values of the form factor V^ Vb \ f with a variation of certain parameters from the lattice NRQCD data on a coarse 
lattice (Set 2 in Table I). Error is statistical only. 


p 

Value 

C(p, U^li^i.oo) 

C(p,V? 1 Vh |f, 04=1 

19) C(p,V£ r,b \F,c 7 = 1.50) C(p, Vo[ ' lb \F, 02=1.25) 

Vii"’ k,=4=1.00 

0.1426(47) 




^2l'| f,c4=1.19 

0.1772(44) 

0.3040 



T/ Tr ?b I 

V 21 \ F,c 7 =1.50 

0.1687(67) 

0.0342 

0.0472 


iz Tr >b 1 

V 21 1^,02=1.25 

0.1769(46) 

0.2979 

0.3352 

0.0467 

■TrTrjb 1 

V 21 \ F,m b =2.59 

0.1939(48) 

0.3070 

0.3479 

0.0508 0.3411 


should be aO(t 2 ) effect. This is roughly consistent. We 
assume a dependence on C 7 as b Cj ~ 2 v 2 b Ci = 0.26 C4 , 
and similarly to the C 4 error above, introduce an additive 
systematic error (correlated across lattice spacings) of 
2 a s v 2 b Ci from not knowing C 7 past tree-level. Just as 
with variations of C 4 , the currents with derivatives are 
insensitive to variations of C 7 and are all consistent within 
our small statistical errors. 


3. Sensitivity and Error from mi,: 

Using the fact that radial splittings are expected to 
be E 21 ~ mbV 2 , by examining (23) we observe that the 
form factor should have a functional dependence on mb 
as V 21 I '’\f = Smi + b mh lm\. Using our data points in 
Table XII, we find a mbl b mb listed in Table XIII. 

Again, we can check consistency within this first or¬ 
der approximation. Comparing the assumed functional 
forms against the equation from which they came (23), 
we find b mb = c\m^b Ci . Thus, using the values of 
6 C4 , b mb we obtain from the lattice data, we find the ratio 
bm h /c\m\b Ci = 0.85(35), consistent with 1.0. 

We allow for a systematic error from the (small) un¬ 
certainty in mistuning the 6 -quark mass estimated from 
[13]. By using the above inverse cubic functional depen¬ 
dence on mb, we find of an error of 3b rnb d mh lm\. Using 
the estimate of b mb in terms of 6 C4 , we find the error as 
3clb c J mh /m b . 


4- Sensitivity and Error from 02 .: 

From our numerical data, it appears as if the form 
factor is not sensitive to a variation in C 2 . We can under¬ 
stand this and use it in our analysis of the errors. In Ap¬ 
pendix C we show how the the leading spin-independent 
perturbative potential from the exchange of a single gluon 
involving the Darwin term at one of the vertices [20] gives 
rise to a correction to the leading order matrix element 
that is 0(a s v 2 ). Using the data in Table XII for how 
V^ Vb \ F varies with C 2 , and using the functional form 


V 2 X Vb |f = CL C2 + C 2 b C2 , we find the values listed in Ta¬ 
ble XIII. 

To test the consistency of this description, by compar¬ 
ing the value b Ci associated with the second term in (23) 
and the second term in (C5) we see b C2 ~ 3v 2 b C 4 /8. Using 
the values in Table XIII gives 3v 2 b Ci /8 = 0.00311(49), 
consistent with b C2 = 0.001(21). Due to the smallness 
of this dependency, we can safely neglect the systematic 
error from not knowing C 2 to two loop order. 


5. Sensitivity and Error from C 3 : 

Since the bottomonium states under study have no or¬ 
bital angular momentum, there is no sensitivity to C 3 
arising from a spin-orbit perturbing potential. This is 
confirmed by the numerical data in Figure 7. We intro¬ 
duce no error from C 3 . 


6. Sensitivity and Error from Four-Quark Operators: 

The four quark operators in NRQCD [13] are contact 
terms between the quark and anti-quark fields arising 
from a 2 processes in relativistic QCD. These can have 
a noticable effect on the hyperfine splitting [16]. Since 
the matrix element in (23) is sensitive to parameters in 
much the same way as the hyperfine splitting, we would 
expect contributions from the four quark operators. In 
Appendix C, we show the effect of the four-quark poten¬ 
tial on the matrix element to first order. 

We introduce a systematic error (correlated across lat¬ 
tice sites) for neglecting these leading order four quark 
operators in our calculation. We estimate this by com¬ 
paring the second term in (23) with the second term in 
(C7) to find an error 27b Ci {d\a s — d 2 a s )/ 16 ir and then 
use the values of d\a s — d 2 a s from [ 20 ] (as corrected per 
[39]). 









17 


TABLE XIII: Values of the functional dependency of V _ 2 i ?i6 | f with parameters from the action using data from Table XII. See 
text for details. Error is statistical only. 


p 

Value 

C(p, a c4 ) C{p,b ci ) C(p, a c7 ) C(p, 6 c7 ) C(p, a c2 ) 

C{p,b c2 ) C(p,a mb ) 

&c4 

0.060(16) 







^c4 

0.083(13) 

-0.974 






a c 7 

0.194(18) 

-0.257 

0.395 





bc7 

-0.017(16) 

0.202 

-0.307 

-0.979 




a c 2 

0.179(24) 

-0.389 

0.510 

0.486 

-0.378 



bc2 

-0.001(21) 

0.366 

-0.459 

-0.403 

0.316 

-0.988 


flmj, 

0.077(34) 

-0.382 

0.487 

0.442 

-0.346 

0.562 

-0.506 

bm b 

2.04(65) 

0.363 

-0.448 

-0.381 

0.300 

-0.512 

0.469 -0.994 


7. Error from Missing Higher Order Operators in the 
NRQCD Action: 

The terms in the action that have not been considered 
are the 0(v 2 ) corrections to the C 2 and C 7 terms. Since 
the coefficient b C2 is already quite small, the v 2 correction 
to this will be negligible within our numerical precision 
and can be neglected. The error from v 2 corrections to 
C 7 is estimated as v 2 b Cj = 2 v 4 b C4 . 


8 . Total Error on V^ Vb \f from Terms in the NRQCD 
Action: 

After performing the final continuum and chiral ex¬ 
trapolation as shown in Section VD, we can obtain a 
breakdown of how each of the uncertainties arising from 
the NRQCD action effects the error in V^ r]b \F as a per¬ 
centage of the error on the total form factor given in 
Table XIV. We find that the errors from the NRQCD ac¬ 
tion contribute to a 10.4% systematic error in V^ Vb \ F as 
a percentage of the total error on the total form factor. 
In order of dominance, the most sizable of these errors is 
a 7.9% error from neglecting the 0(a 2 ) correction in C 4 , 
then a 4.4% error from the statistical error in while 
3.9% comes from neglecting the one-loop correction to 
C 7 . These numbers should be added in quadrature and 
each is a percentage of the total error on the total form 
factor. 

Note that due to the destructive interference between 
the leading order form factor, R 2 i If, and the other 
currents as shown in Section V, the error coming from 
V 21 Vb |f as a percentage of the total error on V^ Vb is larger 
than the errors on V^^If alone. As a result, improve¬ 
ment of the errors coming from the NRQCD action has 
an appreciable effect. 


9. Test of Uncertainties from the NRQCD Action: 

To ensure that we have performed a reasonable esti¬ 
mation of the errors arising from the NRQCD action, we 
have also tuned C 4 against the T(15) — r] b (lS) hyper- 
fine splitting on the coarse lattice denoted set 2 in Table 
I. In a perturbative framework as described above, the 
hyperfine splitting can be pictured as a result of pertur¬ 
bative potentials shifting the unperturbed energies. The 
most sizable of these is the leading order c 2 potential, 
as described in Section VB1, and then the four-quark 
potential, as described in Section VB 6 . In a numerical 
calculation with no four-fermion operators, tuning the 
numerical hyperfine splitting against the experimental 
one would have the effect of absorbing the above four- 
fermion term (among others) into the tuned C 4 . Stated 
more concretely, 

(cf) 2 -7 (c“) 2 d 2 )a s . (25) 

Then, putting (c ! 11 ” 6 ^) 2 into (23) gives exactly the four 
fermion term which we need in (C7). As such, using 
c tuned numer i ca Jiy would include the effect of the four 
fermion operator for this decay automatically. For the 
nonperturbative tuned c error budget, there are no 
C 7 , leading order four-quark, or missing v s operator errors 
as these will be absorded into the value of c •*£ ined and fed 
back into the matrix element calculation automatically. 
However, from (C7) we see there is still an additive sys¬ 
tematic error of 3 u 2 ( 27 / 167 r)a s & C4 from only knowing the 
difference (d\ — d 2 ), and not d\ and d 2 individually. 

The Particle Data Group average for the hyperfine 
splitting is A exp - = 62.3(3.2) MeV [7], while our lat¬ 
tice calculation with C 4 = 1.23 gives A lat = 62.54(46) 
MeV (statistical and scale setting error only). We get a 
value of C 4 Une = 1.230(5)(31) from tuning C 4 against the 
experimental hyperfine splitting, where the first error is 
from the lattice, and the second from experiment. The 
change from the one-loop perturbative value 1.19 to the 
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nonperturbatively tuned 1.230(5)(31) is well-accounted 
for in the error budget (see Sec. VB 8 ) from the statis¬ 
tical error on 5c^ 1 alone, and including the higher order 
corrections to C 4 and the four-quark error is significantly 
over-compensating for this change. 

Rerunning the computation of the form factor with 
C 4 = 1.23, gives a value of V.J t rih = 0.097(14). This in¬ 
cludes all errors, and the only difference from the above 
error budget is that the error in \ f now comes from 
C 4 Une and the error from knowing only the difference 
d 2 — d\. This value is to be compared with the form 
factor from a perturbatively tuned C 4 shown in Section 
VD, i.e., V?™ = 0.089(22). These are entirely consis¬ 
tent, giving evidence that our error budget is a reliable 
estimation of the errors. 

The four-quark operators appear to increase the value 
of the form factor, in a similar way as they do for the 
hyperfine splitting. However, it was found that including 
the four-quark operators in the calculation of the hyper¬ 
fine splitting largely changed the slope of the continuum 
extrapolation but did not shift the final result away from 
the value computed without the four-fermion operators 
included [16]. 

Based on our analysis, we estimate that by tuning C 4 
against the hyperfine splitting on all ensembles and re¬ 
doing the full calculation, one could reduce the error on 
V 21 ’ lb |f to ~ 4%. Also, we estimate that such a calcula¬ 
tion would give an error on the final form factor of ~ 11 % 
(compared against the value given in Table XIV), where 
now the uncertainties in order of dominance are from the 
neglected currents, neglecting the mixing down in uip\ 
and neglecting the one-loop correction to u>wi- 


C. Errors from Missing Higher Order Currents 


Since we are using an effective field theory to study 
this transition, there will be higher order currents which 
we have not included in this study but that contribute to 
the final form factor. The most sizable current which we 
have not included is the D 2 correction to Jwi- Therefore, 
we include a systematic uncertainty (correlated across all 
lattice sites) of v 2 V 2 \ ,b \w\- 


D. Full Error Budget 


After the analysis performed in the previous sections, 
we are now in a position to give a full error budget for 
the form factor V 2 ^ r>b . To compare to experiment, we per¬ 
form a simultaneous lattice spacing and sea quark mass 
extrapolation. We fit results from all ensembles to the 


TABLE XIV: Full error budget for the total form factor V^ Vb 
relevant for the T(2 S) —> 775 (IS')'/ decay from Figure 10. A 
discussion of the uncertainties in V^/ lb \ p is given in Sec. V B 8 . 
The form factor inferred from experimental data in Section II 
is VoT^^lexp = 0.069(14) and has a relative error of 19.74%. 


Error % 

V Tr > b 

V 21 

Systematic 

10.36 

Stats in V^( Vb 

5.48 

Radiative a 2 in wp 

0.83 

Radiative a s in iowi 

4.71 

Radiative a a in u>s 

2.36 

Radiative ct s in cusi 

0.51 

Mixing down in 

3.92 

Missing currents 

7.08 

afm scale 

1.07 

Experimental masses 

0.03 

Priors 

4.18 

Total 

15.81 


form [ .3, 40] 


V( CL , —V]}hys X |^1 

+ ^ ( aA) 2j kj (l + kjxSxm + k j2 (Sx m ) 2 ) 
i= 1,2 


+ 2li6m (l + l 2 (aA) 2 ) 


(26) 


The lattice spacing dependence is set by a scale A = 
500 MeV, Sx m = (ami, — 2.7)/1.5 allows for a mild 
dependence on the effective theory cutoff anib , and 
Sxi = (ami/am s ) — (ami/am s ) p h ys for each ensemble 
with (m;/m s )ph ys = 27.4(1) is taken from lattice QCD 
[41]. We take a Gaussian prior on the leading order a 2 
term to be 0.0(3), as the HISQ action is correct through 
0 (a s a 2 ); a prior of 0 . 0 ( 1 . 0 ) on the higher order a terms; 
a prior of 0.00(3) on allowing for a 3% shift if the light 
quarks were as heavy as the strange; a prior of 0.10(5) 
on Vphys 6 - The extrapolation with all errors is shown in 
Figure 10 and a full error budget is shown in Table XIV. 

By studying the error budget we see that the main 
sources of error are from the systematics in V 2 ^ rib \ f- Here, 
as discussed in Sec. V B 8 , the main sources of uncertainty 
come from the statistical error in C 4 and from not know¬ 
ing the coefficient of a 2 in the expansion of C4. While the 
statistical error on could potentially be reduced from 


The width on this prior is chosen so as to ensure that the fitted 
result is insensitive to the central value. 
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FIG. 10: Fit results for the form factor relevant to the 

T(2S I ) —> ? 7 i,(lS ')7 decay. All errors included. The error bud¬ 
get is shown in Table XIV. 


7% — 10% to 2% — 3% [20] , computation of the two-loop 
coefficient of a 2 would be difficult and lengthy, and un¬ 
likely to be done in the near future. Alternatively, one 
could tune C 4 against the hyperfine splitting on all ensem¬ 
bles, as shown in Section VB9, and the error on V^ Vb 
could be reduced to ~ 11 %. 

After this, the main uncertainty comes from the miss¬ 
ing currents. These could be included with more compu¬ 
tational time if neccessary. While the statistical error on 
each current alone is around 3%, these statistical errors 
do not allow the correlations between the data points in 
the fit to constrain the final result as much as we would 
like, and the final error from statistics in the error budget 
is 5% as a result. Reducing the error from statistics is 
unlikely to have a sizable effect. 

Based on our analysis, we estimate that by including 
the next order of relativistic corrections to the current, 
the mixing down in uip\ and tuning C 4 against the hy¬ 
perfine splitting on all ensembles, an error on V^ nb of 
8 % could be possible (compared against an error of 19% 
on the value inferred from experiment), where the uncer¬ 
tainties in order of dominance would be from the one- 
loop corrections to % 1 and ws and the systematic error 
coming from V^i' lb \ f- 

Our final answer for the form factor is: 

(q 2 = 0) = 0.081(13) (27) 
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FIG. 11: Comparison of our result for the branching frac¬ 
tion (square) with experiment (vertical gray band) and po¬ 
tential model estimates from [1] (crosses). The y-axis labels 
the different references [42-46] and more information about 
these can be found in [ ]. Using the pNRQCD decay rate 
[ ], combined with the experimental total width from the 
PDG average given in Section II, gives a branching fraction 
of 1.9l?;g x 1CT 4 . 


• Previous work only had one lattice spacing. We 
use five ensembles with a fully 0 (a s a 2 ) tadpole- 
improved Liischer-Weisz gluon action with HISQ 
u, d, s and c quarks in the sea, provided by the 
MILC collaboration. These ensembles each have 
~ 1000 configurations and one has physical light 
quark masses. 


Final values for the decay rate and branching fraction are 
given in Section VI. 

VI. DISCUSSION AND CONCLUSIONS 


• We use three relativistic corrections to the leading 
order current as described in Section IV and we 
also test the sensitivities of the form factors from 
all these currents to the parameters in our action 
as shown in Figure 7. 


In this paper we have computed the hindered Ml 
T(2S) —> 7 ^ 5 ( 15)7 decay rate using a lattice NRQCD for¬ 
malism for the 6-quark. We include several improvements 
on earlier exploratory work [5, 6] which are fundamental 
to obtaining an accurate value for this decay rate. The 
key improvements are: 


• We use 0(a s ) correct values for the matching co¬ 
efficients in the NRQCD action. We also take into 
account issues in tuning the 6-quark mass as de¬ 
scribed in Section IIIB. As shown in Figure 7, this 
decay is very sensitive to a subset of these param¬ 
eters. 
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• We calculate the 0(a a ) contribution to the 
matching coefficient of the leading order 
ijj'cr ■ BqedV’/2to(, current which mediates this 
decay, as described in Section IV A. 

• While previous work extracted the matrix element 
by extrapolating/interpolating to the |q| p hys point, 
which only gives the photon on-shell contribution 
q 2 = 0 if the hyperfine splitting is correct, we use 
twisted boundary conditions to extract the form 
factor relevant to this decay at the physical q 2 = 0 
point. 

In Section III E we performed an analysis of the energy 
eigenstates of NRQCD at non-zero momentum. This is 
necessary as the energy eigenstates of a rotationally in¬ 
variant theory, like NRQCD, in an infinite volume con¬ 
tinuum at non-zero momentum are classified by helicity, 
unlike in a Lorentz invariant theory where they are de¬ 
scribed by the standard angular momentum J. This has 
important consequences for a lattice NRQCD calculation 
as additional states appear in the spectrum at non-zero 
momentum (see Figure 1) and one has to be careful to en¬ 
sure that the correct matrix elements are extracted from 
the correlator data. 

In Section V, we show results for the four form fac¬ 
tors from the currents listed in Section IV which when 
renormalised, summed and extrapolated to the contin¬ 
uum limit, can be compared to the form factor inferred 
from experimental data. We found that relativistic cor¬ 
rections to the leading order current gave a negative con¬ 
tribution causing destructive interference, that the power 
counting of the currents deviated from what one would 
naively expect in NRQCD, and that a range of sensitivi¬ 
ties needed to be explained. 

In Section VB, using a simple potential model, we ex¬ 
plained that the matrix element of the leading order cur¬ 
rent was suppressed due to the orthogonality of the radial 
wavefunctions, and this causes the matrix element to be¬ 
come sensitive to a multitude of effects such as relativis¬ 
tic corrections to the leading order current, and certain 
parameters in the NRQCD action that give rise to per¬ 
turbing potentials causing relativistic corrections to the 
wavefunctions, particularly those which effect the hyper¬ 
fine splitting. 

It has been suggested [5, ] that the large changes ex¬ 
perienced in going from an unimproved calculation to an 
improved calculation may mean that it would be benefi¬ 
cial to avoid nonrelativistic approximations. We come to 
a different conclusion and illustrate that although such a 
calculation is intrinsically difficult, NRQCD does indeed 
show that a systematic approach works while also giving 
insight into the process under study. 

After performing the continuum and sea quark mass 
extrapolation, we obtain the form factor V^ 1 b (0)|i at = 
0.081(13), with a full error budget in Table XIV. This 
form factor can be combined with the experimental 


masses used in Section II to produce the decay rate: 

T lat (T(25) -»■ 776(15)7) = 1.72(55) x 10" 2 keV (28) 

which can be compared against the experimental decay 
rate T ex P(T(25) -> 775 ( 15 ) 7 ) = 1.25(49) x 10~ 2 keV [2, 
7]. Using the experimental total width from the PDG 
average given in Section II with our decay rate gives a 
branching fraction of £>(Y(25) —► 776 ( 15 ) 7 ) = 5.4(1.8) x 
10~ 4 which can be compared against the BaBar result of 
3.9(1.5) x 10 -4 [2]. A comparison of our calculation with 
potential model results including relativistic corrections 
[I] is shown in Figure 11. 

Potential model predictions of hindered Ml decay rates 
are known to be particularly difficult to pin down [38] and 
can mischaracterise the experimental data by an order of 
magnitude without relativistic corrections [ 8 ]. Accord¬ 
ing to the Quarkonium Working Group reviews [8, 38], 
sources of uncertainty that contribute to making such de¬ 
cays complicated to calculate include the form of the long 
range potential chosen, and the results depending explic¬ 
itly on the quark mass and the perturbative potential 
chosen. Without relativistic corrections, the branching 
fraction of the Y(25) —► 77 ( 15)7 decay from potential 
model predictions ranges from (0.67 — 11.0) x 10~ 4 [1]. 
Due to the suppression mentioned above, the value of the 
decay rate is very dependent on good knowledge of the 
relativistic corrections [ ]. Including relativistic correc¬ 
tions, potential model predictions for the same branching 
fraction have a wider range (0.05 — 15.0) x 10 -4 , show¬ 
ing indeed that the decay rates may be sensitive to small 
details of the potential [1]. 

The Y(25) —>• ?7b(15)7 decay is sensitive to many of 
the same effects as the hyperfine splitting and an accu¬ 
rate calculation of this decay relies on having the correct 
hyperfine splitting. Given the large range of estimates of 
the hyperfine splitting from potential model predictions 
(46 — 87 MeV [38]), we should not be surprised that the 
potential model estimates for this decay rate also have a 
large range. 

Additionally, radiative transitions between bottomo- 
nium states provide a search for new-physics effects, 
seperate from the weak-sector searches common in the 
literature [ ]. For example, the hyperfine splitting be¬ 

tween the T(15) and ?76(15) has been an important quan¬ 
tity in bottomonium physics, being historically difficult 
for both experimentalists and theorists to predict reli¬ 
ably. Using hindered Ml decays, the BaBar [2, 48] and 
CLEO [ 19] experiments inferred this hyperfine splitting 
to be = 69.3 ± 2.8 MeV [50]. However, in 2012, 

BELLE measured the hb(2P, IP) —» 776 ( 15)7 branching 
fractions (called El decays in the literature), removing 
the dependence on hindered Ml decays and used a signif¬ 
icantly larger sample of events, to yield a hyperfine split¬ 
ting of = 57.9 ± 2.3 MeV [ 1], where A^f' — Ag 3 ) 13 ’ 
has a 3.2cr tension with being zero. 

It has been suggested that the tension of A^* 13 ' and 
theory [16] with A^] 3 ' could, if it persists, indicate a hint 
at new physics [52, 53]. For example, in a multiple-Higgs 
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extension to the standard model, one would speculate 
that the % xp ' seen in experiments is actually an admix¬ 
ture of the true rj b and a CP-odd Higgs boson with mass 
ranging from 9.4 — 10.5 GeV. A relatively light CP-odd 
Higgs scalar can appear in non-minimal supersymmetric 
extensions of the standard model, such as the next-to- 
minimal supersymmetric standard model [53]. In such 
cases, the measured decay rate for T(25) —> 775(15)7 
would likely differ from the Standard Model prediction. 
As stated above, this decay is sensitive in much the 
same way as the hyperfine splitting. To observe a simi¬ 
lar tension between theory and experiment here as that 
existing between Ap* p ' and A^ p ’ would require a 5% 
uncertainty on the form factor (~ 10% on the decay 
rate). The error on the lattice form factor could be re¬ 
duced to ~ 8% (as discussed in Section VD) if more 
precise experimental results became available. Any hint 
of new physics arising from a deviation between the ex¬ 
perimental Y(25) —► 775(15)7 decay rate and theory 
could then be explored more concretely. Additionally, 
the r]b(2S) -7 Y(15)7 decay is an alternative approach 
to studying such effects and a study of this decay rate is 
already underway. 

El radiative decays are more easily computed than hin¬ 
dered Ml decays, and so the El decay rates h b (lP) -7 
775(15)7 and h b {2P) -7 775(15)7 could be calculated 
within this NRQCD framework. Additionally, El cur¬ 
rents can be readily renormalised nonperturbatively. 
Combined with the experimental branching fraction of 
these decays [51], this could give a prediction of the total 
width of the h b (lP) and h b [2P). 
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Appendix A: Classification of Particle States 

Theoretically, particle states living in the Hilbert space 
are classified by unitary irreducible representations (ir- 
reps) of the symmetry group of a theory. We need to con¬ 
sider two symmetry groups here: the Lorentz group and 
the continuous rotational group in three dimensions (the 
symmetry group of NRQCD). The standard procedure to 
build infinite dimensional unitary irreps of these groups 
is via the method of induced representations, where one 
considers finite dimensional unitary irreps of the little 


group and then uses these to build unitary irreps of the 
full group. 

The Poincare group is the symmetry group of a rela¬ 
tivistic quantum field theory, and is given by the semi- 
direct product of the Lorentz group and four transla¬ 
tions. For massive irreps of the Poincare group, the little 
group is 50(3) 7 [5 ]. Thus in a Lorentz invariant the¬ 
ory, massive irreps are defined as \p 2 ; J, M). Note that 
for quarkonia these states are eigenvectors of the charge- 
conjugation operator and parity is also a conserved quan¬ 
tum number, 8 giving the standard \p 2 ; J PC ,M) decom¬ 
position. This description classifies experimental states 
seen to date [7]. 

In a continuum theory that is only rotationally invari¬ 
ant, the analogue of the Poincare group is the semi-direct 
product of the rotational group 50(3) with the three 
translations. For a rotationally invariant theory with zero 
momentum, the little group is also 50(3) and the states 
are classified as |p 2 ; J,M). Thus states in a rotationally 
invariant theory at rest overlap with those in a Lorentz 
invariant theory at rest, where again, parity and charge 
conjugation are good quantum numbers in similar situa¬ 
tions. Given that at nonzero momentum in a rotationally 
invariant theory we cannot perform a Lorentzian boost 
to the rest frame, the little group at nonzero momentum 
is now different to the zero momentum little group. The 
little group is now 50(2) 9 [■ ]. In this case, the uni¬ 
tary irreps are classified by |p 2 ; A), where A is an eigen¬ 
value of the hclicity operator A = p ■ J/E. The helicity 
A = Ao will get contributions from all J with Ao < J. 
This can have important consequences for the extracted 
energy spectrum in NRQCD, c.f., Figure 1 and 2, and is 
fundamentally different from the Lorentzian theory. 

At zero momentum, the operators 77 s and 7* that we 
use in this calculation overlap onto 0 p and 1 states in 
a rotationally invariant continuum theory [26]. We now 
need to find which helicity eigenstates these operators 
overlap with at nonzero momentum. The authors of [26] 
illustrate how to construct helicity operators via 

© J ’ A (P) = E ^mx(R)0 JM (p) (Al) 

M 

where T>^ IX (R) is a Wigner-P matrix, R is the active 
transformation which rotates (0,0, |p|) to p, 0 J,A (p) is 
a helicity operator with helicity A in an infinite volume 


7 At nonzero momentum we can perform a Lorentz boost back to 
the rest frame, ensuring the little group is the same for zero and 
nonzero three-momentum 

8 At nonzero momentum, these states are not eigenvectors of the 
parity operator, but are eigenstates of the n operator defined in 
the text, which conserves parity. 

9 The construction of the irreps for a rotationally invariant theory 
at nonzero momentum is similar to a massless representation in 
a Lorentz invariant theory. 
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continuum, e.g., 

(0|O J > A (p)|p; J', A') = Z^ J '^S xy (A2) 

and we refer the reader to Ref. [26] for further de¬ 
tails. For quarkonium, the possibile values of A = 
{0+,0“, |1|, |2|,...}, where the +/— on the A = 0 rep¬ 
resent the II symmetry with eigenvalue fj = P(—1) J [26]. 
Using the fact that the Wigner-2? matrices with J = 0 
are 5 X m, the O 7 , O' 1 ' bilinear operators which we use 
in this calculation give rise to the helicity operators at 
nonzero momentum 

O J= 0’ A=0 '(p) = eU 5 ( p) 

O J=i, A =o+ (p) = J2v J m = x %(R)CP M (p) 

M 

0 j=i,a=|*I( p ) = £ V j ^1 w {R)0^ M (p). (A3) 

M 

As can be seen, O^ ( p) is a helicity operator which cre¬ 
ates a A = 0“ state, but O 7 (p) creates an admixture of 
A = 0 + , |1| states. 

The question now is: how do we identify which J PC 
contributes to each A, and how do we parameterise the 
amplitudes? By noticing that the helicity A = J z when 
the momentum is projected onto the z-axis, all states 
with J > A will have a J. large enough to give a contri¬ 
bution to this helicity state (see Table III). 

We also want to know how to quantify the amplitudes. 
In a rotationally invariant theory, the invariant quantities 
are Sij and £ijk- For a J p state, we also have the momen¬ 
tum pj and the symmetric polarisation tensor 
We can use these to parameterise the amplitudes rele¬ 
vant for a rotationally invariant theory. For the operator 
O 7 , Table XI in [26] has the possible decompositions and 
we reproduce the parameterisations for the O 7 operator 
which are important for our calculation 

(0\O^(p)\n0- + (p)) = Z n (A4) 

(0|O 7 ' J (p)|nl ++ (e,p)) = Z' n e i p i /m nl ++ 

(O|0 7 ' J (p)|n2 f (e,p)} = Z^eu +Zle i:i pip j /m 2 n2 - + . 

where n is the radial label. Using the overlap for the 
1 ++ from (A4) to parameterise the continuum two-point 
correlator with nonzero momentum, one finds that the 
amplitudes from our fit with local smearing should be 
suppressed by |p|/mi++ relative to states which overlap 
with the operator at zero momentum. For the momen¬ 
tum that we use in our calculation, this factor is around 
7%, and we observe that in our correlator data, the am¬ 
plitudes for the states which do not overlap at zero mo¬ 
mentum (and for which we get a signal) such as the 1 ++ , 
are suppressed by this factor while the other amplitudes 
are 0(1). We observe that as the momentum increases, 
so does the value of the amplitude at fixed lattice spacing. 

Additionally, the symmetry group giving rise to the 
invariants which classify states, e.g., the little group, is 


broken by a finite volume lattice to a reduced symmetry 
group [55]. At zero momentum with a cubic lattice, this 
reduced symmetry group for quarkonia is the octahedral 
group, Oh- States are now classified in terms of irreps 
of Oh, denoted A PC , where [56] shows how to subduce 
operators with continuum spin J PC to operators with 
definite A PC on the lattice. As mentioned above, in an 
infinite-volume continuum theory, the O 7 ( O 7 ) opera¬ 
tor overlaps only with J PC = 0 ^(l ) at rest, but this 

operator falls into the A]~ + (Tj ) irrep of Oh on the lat¬ 
tice, where mixing with the J pc = 4 h (3 ) state (and 

higher spins) is possible. However we do not see this mix¬ 
ing: rotational symmetry breaking is found to be weakly 
broken with a fine lattice and with a rotationally invari¬ 
ant smearing for a particular lattice setup [56], where 
the spectrum and overlaps were compatible with an effec¬ 
tive restoration of rotational symmetry. For this reason, 
we choose to use a rotationally invariant smearing, an 
isotropic lattice and have discretisation improvements in 
our action. Secondly, the masses of the additional states 
are too large to be seen in the first few energy levels which 
we are interested in. As such, they will only potentially 
contribute as additional discretisation effects in the low¬ 
est energy modes. Indeed, studies of the spectrum from 
NRQCD by the HPQCD collaboration indicate this to 
be the case (see Appendix C of [13]). 

For the nonzero momentum case, the reduced little 
group actually depends on the type of momenta. This 
is due to the fact that a general integer-valued momen¬ 
tum on the lattice cannot be rotated into the z-axis like 
in an infinite volume continuum, 10 e.g. there is no oc¬ 
tahedral transformation which rotates (0,1,1) to the z- 
axis. We use an isotropic momentum (rather than an 
on-axis momentum) as it has been shown to break rota¬ 
tional invariance less and lead to smaller discretisation ef¬ 
fects [13]. For our isotropic momentum, the reduced little 
group is Dic 3 [26] . The operator O^ (O 7 ’ ) falls into the 
A 2 {A\ and E 2 ) irrep of Dic 3 , where mixing with A = 3 
(3 and 2) states is possible. For O 7 , this gives rise to 
potential mixing from 3 ±+ , 4 ±+ states (and higher spin). 
As in the zero-momentum case, this mixing due to the 
lattice was found to be negligible with a fine lattice and 
a rotationally invariant smearing for a particular setup 
[26]. These states should be of higher energy than the 
first few states in our spectrum, and we see no evidence 
of them in our low lying spectrum. For the O' y opera¬ 
tor, there can be mixing with A = 2 (2 < J with J z = 2 
states) which is not important for our analysis. 

There is an important distinction to be understood 
from using a rotationally invariant formalism for the 
quark versus a Lorentz-invariant one. If each of these 


10 With twisted boundary conditions, the momenta are still discre- 
tised but just shifted by an arbitrary value. As such, the little 
group of momentum with a twist is the same as the little group 
of momentum without a twist. 
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formalisms is discretised, then at fixed nonzero momen¬ 
tum, the discretised version of the Lorentz-invariant the¬ 
ory might be broken to a rotationally invariant theory, 
e.g., by using an anisotropic lattice spacing in the time 
direction. As such, as the infinite volume continuum 
limit is taken, any overlap onto J PC as a result of helic- 
ity eigenstates (such as the 1 ++ from the O 7 operator) 
would disappear [57]. However, in a rotationally invari¬ 
ant theory like NRQCD, as the lattice spacing is taken 
to zero, these overlaps are still present as they are an 
infinite volume continuum effect. This is why we find a 
similar signal across all lattice spacings for these states 
in NRQCD. 


Appendix B: Twisted Correlators with Derivative 
Operators 


For clarity, we will describe the construction of the 
twisted correlators with derivative operators in this sec¬ 
tion. To gain access to arbitrary momenta on the lattice, 
one can define a quark field [22, 23] that satisfies 9 BC 
via x + etL ) = e l2 ' nXi 'ip 0 (x), where 9i = 2-nxi/L. Now 
the available momentum space is A = {k = p + 0\ki = 
2n(rii + Xi)/L, where n, € Z}. Notice that the available 
momentum space has an arbitary shifted value 6 that 
we can choose to give the physical point q 2 = 0. One 
now builds interpolating operators from these 9 BC fields 
as 0(x;9 2 9i) = %j) 02 {x)T‘ip 01 {x), which gives rise to the 
two-point correlator 


C 2pt (0i - e 2 + p, t) = J2 e-^- e2+p >- x 


Tr 


(T^ 2 (0,0|x, t))(T f S 01 (x, i|0, 0)) (Bl) 


where 5 e (0,0|x, f) is a quark propagator found by in¬ 
verting the Dirac matrix, D e (x,y), defined via S[il> 0 } = 
J2 X y w (%)D e (x, y)ip 0 (y). As a consequence of ip 6 sat¬ 
isfying 0BC, the Dirac matrix D 0 {x,y) also satisfies the 
same boundary conditions. This is an inconvenience as 
typical inverters are built with PBC. However, it is pos¬ 
sible to use a trick in order to use the PBC invertors yet 
still get access to the 0BC correlator data in (Bl). 

To do this, one notices that a second quark field, de¬ 
fined via the scaling ip 0 {x) = e~ 2 ' Kl0 ' Xl,L 'ijj 0 (x), satisfies 
PBC yet still includes information on the twist. Since 


S e (x\y) = e ie -(*-y ) S 0 (x\y) (B2) 


S 0 (x\y) is a quark propagator found by invert¬ 
ing the Dirac matrix, D 0 (x,y), where D 0 (x,y) = 
e- w ' x D 9 (x,y)e l6 ' y . D e (x\y) satisfies PBC by construc¬ 
tion and the two exponentials only alter the derivative 
in the Dirac action and can be implemented by scaling 
the gluonic fields (before inverting) as C7 /i ( x) —> U 0 (x) = 
e i2n / L0 » Utl (x) with (9 m = (0, 6 ) [22], 

The final step is to rewrite the twisted correlator in 
(Bl) in terms of the propagator we actually compute us¬ 
ing (B2) 


C 2pt (0i - e 2 + p,i) = Y, e-' (ei - 02 + p) - x 

X 

x Tr [(rje _i6,2 ' x 5' 6 ' 2 (0, 0|x, t)) ((Tfe^^S 01 (x, i|0,0))] . 

(B3) 


If T = V, then 


C 2 pt(0i -02+P,i) =^e- ip - x '&[(e ie2 - x V fe e- i02 - x 5 02 (O,O|x,t)) {e~ i0 ^\7 k e i6 ^S 0 ' (x,i|0,0)) ] . (B4) 

X 


This can be implemented in the same way as the twist 
in the Dirac invertor, by using U 0 ( x ) in the construction 
of the covariant derivative operator. This “changing the 
derivatives” issue does not occur in our two-point corre¬ 
lators, but does occur in the (more complicated) three 


point correlators with currents Jwi, Js, Jsi from (16). 
To give an explicit example of the three point correla¬ 
tor using the current Jwi, by keeping the initial state at 
rest, and twisting only one propagator in the final state 
with 9f, we have 


^3pt(Pf = Pf + Q t , q e = q - 0 t ; t, T) = -i ]T e _ip f x Tr [s 0 f (x, T|y, t) 

x,y 

|d 2 , (cr X q e Q n e- j(q - 0f) - y }s , (y,t|O,O)") CT m S(0,0|x,T) 


8m 3 b 


(B5) 
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where we can clearly see that D 2 does not commute with 
e -ie f -y, but not all derivatives are twisted due to the 
commutation. Since there are no derivatives in the Jp 
current, the Of terms cancel and this issue is avoided. 
Smearing the twisted fields leads to a similar issue as 
presented above with the derivative, and so we do not 
smear the twisted fields. Analogous complications arise 
when using point-split operators with twisted momentum 
in staggered quark formalisms [58] . If done correctly, and 
any smearings are applied appropriately, the correlator 
data from using dBC and PBC should agree on a configu¬ 
ration basis to machine precision (if the total momentum 
is identical for all states). 

Appendix C: Error Analysis Using a Simple 
Potential Model 

First, we want to find the sensitivity of the matrix ele¬ 
ment to c\ using a potential from the exchange of a single 
gluon between two vertices involving the chromomagen- 
tic operator [20]. We find (assuming the wavefunctions 
at the origin for the ij b and T are the same) 

v™ = °)vao) 

v L = 2 ^r n mm(o). (ci) 

Putting this back into (22) with the Jp current yields: 


Then substituting this back into (22) we find: 

^{ Vb (is)\j F \r( 2 S)}^ = 

W( Vb (is)\j F \r(2S))W 

- gfj (e ^T" 10) <0 >(*(mS)|*|T(2S))(»> 

k \ m^l m 

+ rjomo) m Mls)l j FlT(nS)r >) m 

En2 J 

= (0) (vb(is)\J F \r(2S))W 
- ^^^(°)V>2(°)( { 0 ) (r, b (2S)\J F |T(2S)><°> 

- (0) (?7 b (l^)kF|T(n5))( 0 )+O(u 2 ^ . (C5) 

Using (20), we see the leading order terms in the sec¬ 
ond piece of (C5) cancel and we are left with 0(a s v 2 ) 
corrections to the unperturbed matrix element. 

The four quark potential is (assuming the wavefunc¬ 
tions at the origin of the two states are the same) [20] 

Od, O' 2 4 

Kt - °) 

VZn = ^^^(0)^(0) • (C6) 


«< % (1S)|J F |T(2S)}« = 

^( Vb (lS)\J F \T(2S))^ + 


Putting this into (22) and performing an identical anal¬ 
ysis as done above gives 


E (»(, l(m s)|j P |T(2s))i») 

^ \ m^l m 

- ]T 2 V’n(0)V’2(0) (0) <??fc(15) | Jp | T (^))(0)^ . (C2) 

/ r> ^n 2 I 

In getting to (C2) we have used the fact that E^ m = E% b m 
as the unperturbated Hamiltonian has no spin terms. We 
have neglected the T (pS) -a rj b (lS) transitions for p > 
2 in the sum due to the fact that the radial overlap, 
(20), is suppressed by at least 0(v 2 ). In fact, they will 
be suppressed more due to the radial difference getting 
larger and the wavefunction at the origin getting smaller 
for higher radial excitations. Eqn. (23) can be found 
straightforwardly by factoring the spin part of the matrix 
element from the radial part, i.e., using (20). 

If we now consider a potential from the exchange of 
a single gluon involving the Darwin term at one of the 
vertices, we find [20] 

2 

v™=v? m = -^Lr n ( 0hM0). (C3) 


^( Vb (lS)\JF\T( 2 S))^ = 

^(r, b (lS)\J F \T(2S))^ 

9cha 2 4 ^ Vh*(0)Va0) (o) (rj b (mS) | |T(26 1 ))(°) 


3771? 


9d20t 2 


2 E n 2 

u 2 

= ^(r ]b (lS)\J F \T(2S))^ 
, 9 4 ^(0)^(0) 


^ r n (v)Mo) io) {vb{ls) \jF\ T{nS)) m 


2 3 ml 


E 2 


S, 


if 


(d 2 a 2 s - dia 2 s 


0((2d 2 a 2 s -d ia 2 s )v 2 )). (C7) 


The error in the last line was introduced by expanding 
out the radial overlap (20) and noting that the two ma¬ 
trix elements do not have to be identical to first order in 
\q-y\ 2 . Even if we did include the four fermion operators 
in the calculation, since only the combination d\ — d 2 is 
currently known perturbatively, and not d\ and d 2 indi¬ 
vidually, we would still need to introduce the 0 {v 2 ) error 
in our calculation. 
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